JKR Model Implementation

See this page for the documentation of this contact model.

• contactmodeljkr.h
• contactmodeljkr.cpp

contactmodeljkr.h

// contactmodeljkr.h
#include "contactmodel/src/contactmodelmechanical.h"

#include "module/interface/itablelist.h"
#include "module/interface/itable.h"

#ifdef JKR_LIB
#  define JKR_EXPORT EXPORT_TAG
#elif defined(NO_MODEL_IMPORT)
#  define JKR_EXPORT
#else
#  define JKR_EXPORT IMPORT_TAG
#endif

namespace cmodelsxd {
	using namespace itasca;

	class ContactModelJKR : public ContactModelMechanical {
	public:
		// Constructor: Set default values for contact model properties.
		JKR_EXPORT ContactModelJKR();
		// Destructor, called when contact is deleted: free allocated memory, etc.
		JKR_EXPORT virtual ~ContactModelJKR();
		// Contact model name (used as keyword for commands and FISH).
		virtual QString  getName() const { return "jkr"; }
		// The index provides a quick way to determine the type of contact model.
		// Each type of contact model in PFC must have a unique index; this is assigned
		// by PFC when the contact model is loaded. This index should be set to -1
		virtual void     setIndex(int i) { index_ = i; }
		virtual int      getIndex() const { return index_; }
		// Contact model version number (e.g., MyModel05_1). The version number can be
		// accessed during the save-restore operation (within the archive method,
		// testing {stream.getRestoreVersion() == getMinorVersion()} to allow for 
		// future modifications to the contact model data structure.
		virtual uint     getMinorVersion() const;
		// Copy the state information to a newly created contact model.
		// Provide access to state information, for use by copy method.
		virtual void     copy(const ContactModel *c);
		// Provide save-restore capability for the state information.
		virtual void     archive(ArchiveStream &);
		// Enumerator for the properties.
		enum PropertyKeys {
			kwShear = 1
			, kwPoiss
			, kwFric
			, kwjkrF
			, kwjkrS
			, kwDpNRatio
			, kwDpSRatio
			, kwDpMode
			, kwDpF
			, kwResFric
			, kwResMoment
			, kwResS
			, kwKsFac				
			, kwSurfAdh
			, kwActiveMode
			, kwA0
			, kwPullOff
			, kwTearOff
		};
		// Contact model property names in a comma separated list. The order corresponds with
		// the order of the PropertyKeys enumerator above. One can visualize any of these 
		// properties in PFC automatically. 
		virtual QString  getProperties() const {
			return "jkr_shear"
				",jkr_poiss"
				",fric"
				",jkr_force"
				",jkr_slip"
				",dp_nratio"
				",dp_sratio"
				",dp_mode"
				",dp_force"
				",rr_fric"
				",rr_moment"
				",rr_slip"
				",ks_fac"
				",surf_adh"
				",active_mode"
				",a0"
				",pull_off_f"
				",tear_off_d";
		}
		// Enumerator for the energies.
		enum EnergyKeys {
			kwEStrain = 1
			, kwERRStrain
			, kwESlip
			, kwERRSlip
			, kwEDashpot
		};
		// Contact model energy names in a comma separated list. The order corresponds with
		// the order of the EnergyKeys enumerator above. 
		virtual QString  getEnergies() const {
			return "energy-strain-jkr"
				",energy-rrstrain"
				",energy-slip"
				",energy-rrslip"
				",energy-dashpot";
		}
		// Returns the value of the energy (base 1 - getEnergy(1) returns the estrain energy).
		virtual double   getEnergy(uint i) const;
		// Returns whether or not each energy is accumulated (base 1 - getEnergyAccumulate(1) 
		// returns whether or not the estrain energy is accumulated which is false).
		virtual bool     getEnergyAccumulate(uint i) const;
		// Set an energy value (base 1 - setEnergy(1) sets the estrain energy).
		virtual void     setEnergy(uint i, const double &d); // Base 1
		// Activate the energy. This is only called if the energy tracking is enabled. 
		virtual void     activateEnergy() { if (energies_) return; energies_ = NEWC(Energies()); }
		// Returns whether or not the energy tracking has been enabled for this contact.
		virtual bool     getEnergyActivated() const { return (energies_ != 0); }

		// Enumerator for contact model related FISH callback events. 
		enum FishCallEvents {
			fActivated = 0
			, fSlipChange
		};
		// Contact model FISH callback event names in a comma separated list. The order corresponds with
		// the order of the FishCallEvents enumerator above. 
		virtual QString  getFishCallEvents() const {
			return
				"contact_activated"
				",slip_change";
		}

		// Return the specified contact model property.
		virtual QVariant getProperty(uint i, const IContact *) const;
		// The return value denotes whether or not the property corresponds to the global
		// or local coordinate system (TRUE: global system, FALSE: local system). The
		// local system is the contact-plane system (nst) defined as follows.
		// If a vector V is expressed in the local system as (Vn, Vs, Vt), then V is
		// expressed in the global system as {Vn*nc + Vs*sc + Vt*tc} where where nc, sc
		// and tc are unit vectors in directions of the nst axes.
		// This is used when rendering contact model properties that are vectors.
		virtual bool     getPropertyGlobal(uint i) const;
		// Set the specified contact model property, ensuring that it is of the correct type
		// and within the correct range --- if not, then throw an exception.
		// The return value denotes whether or not the update has affected the timestep
		// computation (by having modified the translational or rotational tangent stiffnesses).
		// If true is returned, then the timestep will be recomputed.
		virtual bool     setProperty(uint i, const QVariant &v, IContact *c);
		// The return value denotes whether or not the property is read-only
		// (TRUE: read-only, FALSE: read-write).
		virtual bool     getPropertyReadOnly(uint i) const;

		// The return value denotes whether or not the property is inheritable
		// (TRUE: inheritable, FALSE: not inheritable). Inheritance is provided by
		// the endPropertyUpdated method.
		virtual bool     supportsInheritance(uint i) const;
		// Return whether or not inheritance is enabled for the specified property.
		virtual bool     getInheritance(uint i) const { assert(i < 32); quint32 mask = to<quint32>(1 << i);  return (inheritanceField_ & mask) ? true : false; }
		// Set the inheritance flag for the specified property.
		virtual void     setInheritance(uint i, bool b) { assert(i < 32); quint32 mask = to<quint32>(1 << i);  if (b) inheritanceField_ |= mask;  else inheritanceField_ &= ~mask; }

		// Prepare for entry into ForceDispLaw. The validate function is called when:
		// (1) the contact is created, (2) a property of the contact that returns a true via
		// the setProperty method has been modified and (3) when a set of cycles is executed
		// via the {cycle N} command.
		// Return value indicates contact activity (TRUE: active, FALSE: inactive).
		virtual bool    validate(ContactModelMechanicalState *state, const double &timestep);
		// The endPropertyUpdated method is called whenever a surface property (with a name
		// that matches an inheritable contact model property name) of one of the contacting
		// pieces is modified. This allows the contact model to update its associated
		// properties. The return value denotes whether or not the update has affected
		// the time step computation (by having modified the translational or rotational
		// tangent stiffnesses). If true is returned, then the time step will be recomputed.  
		virtual bool    endPropertyUpdated(const QString &name, const IContactMechanical *c);
		// The forceDisplacementLaw function is called during each cycle. Given the relative
		// motion of the two contacting pieces (via
		//   state->relativeTranslationalIncrement_ (Ddn, Ddss, Ddst)
		//   state->relativeAngularIncrement_       (Dtt, Dtbs, Dtbt)
		//     Ddn  : relative normal-displacement increment, Ddn > 0 is opening
		//     Ddss : relative  shear-displacement increment (s-axis component)
		//     Ddst : relative  shear-displacement increment (t-axis component)
		//     Dtt  : relative twist-rotation increment
		//     Dtbs : relative  bend-rotation increment (s-axis component)
		//     Dtbt : relative  bend-rotation increment (t-axis component)
		//       The relative displacement and rotation increments:
		//         Dd = Ddn*nc + Ddss*sc + Ddst*tc
		//         Dt = Dtt*nc + Dtbs*sc + Dtbt*tc
		//       where nc, sc and tc are unit vectors in direc. of the nst axes, respectively.
		//       [see {Table 1: Contact State Variables} in PFC Model Components:
		//       Contacts and Contact Models: Contact Resolution]
		// ) and the contact properties, this function must update the contact force and
		// moment.
		//   The force_ is acting on piece 2, and is expressed in the local coordinate system
		//   (defined in getPropertyGlobal) such that the first component positive denotes
		//   compression. If we define the moment acting on piece 2 by Mc, and Mc is expressed
		//   in the local coordinate system (defined in getPropertyGlobal), then we must use the getMechanicalContact()->updateResultingTorquesLocal(...) method to 
		// get the total moment. 
		// The return value indicates the contact activity status (TRUE: active, FALSE:
		// inactive) during the next cycle.
		// Additional information:
		//   * If state->activated() is true, then the contact has just become active (it was
		//     inactive during the previous time step).
		//   * Fully elastic behavior is enforced during the SOLVE ELASTIC command by having
		//     the forceDispLaw handle the case of {state->canFail_ == true}.
		virtual bool    forceDisplacementLaw(ContactModelMechanicalState *state, const double &timestep);
		// The getEffectiveXStiffness functions return the translational and rotational
		// tangent stiffnesses used to compute a stable time step. When a contact is sliding,
		// the translational tangent shear stiffness is zero (but this stiffness reduction
		// is typically ignored when computing a stable time step). If the contact model
		// includes a dashpot, then the translational stiffnesses must be increased (see
		// Potyondy (2009)).
		//   [Potyondy, D. 'Stiffness Matrix at a Contact Between Two Clumps,' Itasca
		//   Consulting Group, Inc., Minneapolis, MN, Technical Memorandum ICG6863-L,
		//   December 7, 2009.]
		virtual DVect2  getEffectiveTranslationalStiffness() const { return effectiveTranslationalStiffness_; }
		virtual DAVect  getEffectiveRotationalStiffness() const { return effectiveRotationalStiffness_; }

		// Return a new instance of the contact model. This is used in the CMAT
		// when a new contact is created. 
		virtual ContactModelJKR *clone() const { return NEWC(ContactModelJKR()); }
		// The getActivityDistance function is called by the contact-resolution logic when
		// the CMAT is modified. Return value is the activity distance used by the
		// checkActivity function.
		virtual double              getActivityDistance() const { return (rgap_ + distance_active_); }		// initially 'distance_active_' is zero so that the contat becomes active when the particles touch (zerp overlap)
																											// when the contact becomes active,it is set so that the contact only becomes inactive at this tear_off_ distance
		// The isOKToDelete function is called by the contact-resolution logic when...
		// Return value indicates whether or not the contact may be deleted.
		// If TRUE, then the contact may be deleted when it is inactive.
		// If FALSE, then the contact may not be deleted (under any condition).
		virtual bool                isOKToDelete() const { return !isBonded(); }
		// Zero the forces and moments stored in the contact model. This function is called
		// when the contact becomes inactive.
		virtual void                resetForcesAndMoments() {
			jkr_F(DVect(0.0)); dp_F(DVect(0.0));
			res_M(DAVect(0.0));																			
			if (energies_) {
				//energies_->estrain_ = 0.0;
				energies_->errstrain_ = 0.0;
			}
			distance_active_ = 0.0;																			// reset 
		}
		virtual void     setForce(const DVect &v, IContact *c);
		virtual void     setArea(const double&) { throw Exception("The setArea method cannot be used with the JKR contact model."); }
		virtual double   getArea() const { return 0.0; }
		// The checkActivity function is called by the contact-resolution logic when...
		// Return value indicates contact activity (TRUE: active, FALSE: inactive).
		// A contact with the arrlinear model is active if the surface gap is less than
		// or equal to the attraction range (a_d0_).
		virtual bool     checkActivity(const double &gap) { return  gap <= (rgap_ + distance_active_); }	// initially 'distance_active_' is zero so that the contat becomes active when the particles touch (zerp overlap)
																											// when the contact becomes active,it is set so that the contact only becomes inactive at this tear_off_ distance
		// Returns the sliding state (FALSE is returned if not implemented).
		virtual bool     isSliding() const { return jkr_S_; }
		// Returns the bonding state (FALSE is returned if not implemented).
		virtual bool     isBonded() const { return false; }

		// Both of these methods are called only for contacts with facets where the wall 
		// resolution scheme is set the full. In such cases one might wish to propagate 
		// contact state information (e.g., shear force) from one active contact to another. 
		// See the Faceted Wall section in the documentation. 

		// Return the total force that the contact model holds.
		virtual DVect    getForce(const IContactMechanical *) const;

		// Return the total moment on 1 that the contact model holds
		virtual DAVect   getMomentOn1(const IContactMechanical *) const;

		// Return the total moment on 1 that the contact model holds
		virtual DAVect   getMomentOn2(const IContactMechanical *) const;

		// Methods to get and set properties. 
		const double & shear() const { return shear_; }
		void           shear(const double &d) { shear_ = d; }
		const double & poiss() const { return poiss_; }
		void           poiss(const double &d) { poiss_ = d; }
		const double & fric() const { return fric_; }
		void           fric(const double &d) { fric_ = d; }
		const DVect &  jkr_F() const { return jkr_F_; }
		void           jkr_F(const DVect &f) { jkr_F_ = f; }
		bool           jkr_S() const { return jkr_S_; }
		void           jkr_S(bool b) { jkr_S_ = b; }
		const double & rgap() const { return rgap_; }
		void           rgap(const double &d) { rgap_ = d; }
		
		const double & kn() const { return kn_; }
		void           kn(const double &d) { kn_ = d; }
		const double & ks() const { return ks_; }
		void           ks(const double &d) { ks_ = d; }
		const double & ks_fac() const { return ks_fac_; }
		void           ks_fac(const double &d) { ks_fac_ = d; }
		const double & surf_adh() const { return surf_adh_; }
		void           surf_adh(const double &d) { surf_adh_ = d; }
		const int    & active_mode() const { return active_mode_; }
		void           active_mode(const int &d) { active_mode_ = d; }
		const double & distance_active() const { return distance_active_; }
		void           distance_active(const double &d) { distance_active_ = d; }
		const double & a0() const { return a0_; }
		void           a0(const double &d) { a0_ = d; }
		const double & pull_off() const { return pull_off_f_; }
		void           pull_off(const double &d) { pull_off_f_ = d; }
		const double & tear_off() const { return tear_off_d_; }
		void           tear_off(const double &d) { tear_off_d_ = d; }
		const double & r_hertz() const { return r_hertz_; }
		void           r_hertz(const double &d) { r_hertz_ = d; }
		const double & c1() const { return c1_; }
		void           c1(const double &d) { c1_ = d; }
				
		bool     hasDamping() const { return dpProps_ ? true : false; }
		double   dp_nratio() const { return (hasDamping() ? (dpProps_->dp_nratio_) : 0.0); }
		void     dp_nratio(const double &d) { if (!hasDamping()) return; dpProps_->dp_nratio_ = d; }
		double   dp_sratio() const { return hasDamping() ? dpProps_->dp_sratio_ : 0.0; }
		void     dp_sratio(const double &d) { if (!hasDamping()) return; dpProps_->dp_sratio_ = d; }
		int      dp_mode() const { return hasDamping() ? dpProps_->dp_mode_ : -1; }
		void     dp_mode(int i) { if (!hasDamping()) return; dpProps_->dp_mode_ = i; }
		DVect    dp_F() const { return hasDamping() ? dpProps_->dp_F_ : DVect(0.0); }
		void     dp_F(const DVect &f) { if (!hasDamping()) return; dpProps_->dp_F_ = f; }

		bool    hasEnergies() const { return energies_ ? true : false; }
		double  estrain() const { return hasEnergies() ? energies_->estrain_ : 0.0; }
		void    estrain(const double &d) { if (!hasEnergies()) return; energies_->estrain_ = d; }
		double  errstrain() const { return hasEnergies() ? energies_->errstrain_ : 0.0; }
		void    errstrain(const double &d) { if (!hasEnergies()) return; energies_->errstrain_ = d; }
		double  eslip() const { return hasEnergies() ? energies_->eslip_ : 0.0; }
		void    eslip(const double &d) { if (!hasEnergies()) return; energies_->eslip_ = d; }
		double  errslip() const { return hasEnergies() ? energies_->errslip_ : 0.0; }
		void    errslip(const double &d) { if (!hasEnergies()) return; energies_->errslip_ = d; }
		double  edashpot() const { return hasEnergies() ? energies_->edashpot_ : 0.0; }
		void    edashpot(const double &d) { if (!hasEnergies()) return; energies_->edashpot_ = d; }
		
		uint inheritanceField() const { return inheritanceField_; }
		void inheritanceField(uint i) { inheritanceField_ = i; }

		const DVect2 & effectiveTranslationalStiffness()  const { return effectiveTranslationalStiffness_; }
		void           effectiveTranslationalStiffness(const DVect2 &v) { effectiveTranslationalStiffness_ = v; }
		const DAVect & effectiveRotationalStiffness()  const { return effectiveRotationalStiffness_; }
		void           effectiveRotationalStiffness(const DAVect &v) { effectiveRotationalStiffness_ = v; }

		// Rolling resistance methods
		const double & res_fric() const { return res_fric_; }
		void           res_fric(const double &d) { res_fric_ = d; }
		const DAVect & res_M() const { return res_M_; }
		void           res_M(const DAVect &f) { res_M_ = f; }
		bool           res_S() const { return res_S_; }
		void           res_S(bool b) { res_S_ = b; }
		const double & kr() const { return kr_; }
		void           kr(const double &d) { kr_ = d; }
		const double & fr() const { return fr_; }
		void           fr(const double &d) { fr_ = d; }
		const double & rbar_square() const { return rbar_square_; }
		void           rbar_square(const double &d) { rbar_square_ = d; }
		//
				
	private:
		// Index - used internally by PFC. Should be set to -1 in the cpp file. 
		static int index_;

		// Structure to store the energies. 
		struct Energies {
			Energies() : estrain_(0.0), errstrain_(0.0), eslip_(0.0), errslip_(0.0), edashpot_(0.0) {}
			double estrain_;   // elastic energy stored in linear group 
			double errstrain_; // elastic energy stored in rolling resistance group
			double eslip_;     // work dissipated by friction 
			double errslip_;   // work dissipated by rolling resistance friction 
			double edashpot_;  // work dissipated by dashpots
		};

		// Structure to store dashpot quantities. 
		struct dpProps {
			dpProps() : dp_nratio_(0.0), dp_sratio_(0.0), dp_mode_(0), dp_F_(DVect(0.0)) {}
			double dp_nratio_;		// normal viscous critical damping ratio
			double dp_sratio_;		// shear  viscous critical damping ratio
			int    dp_mode_;		// for viscous mode (0-1) 0 = no limit or cutoff;  1 = shear cut-off when sliding
			DVect  dp_F_;			// Force in the dashpots
		};

		bool   updateEndStiffCoef(const IContactMechanical *con);
		bool   updateEndFric(const IContactMechanical *con);
		bool   updateEndResFric(const IContactMechanical *con);
		bool   updateEndSurfAdh(const IContactMechanical *con);
		void   updateStiffCoef(const IContactMechanical *con);
		void   updateEffectiveStiffness(ContactModelMechanicalState *state);
		void   setDampCoefficients(const double &mass, double *kn_tangent, double *ks_tangent, double *vcn, double *vcs);

		bool SetJKRTable(ContactModelMechanicalState *state);
		DVect2 InterpolateJKRTable_Method1(double overlap_ratio);
		DVect2 InterpolateJKRTable_Method2(double overlap_ratio, double trans_n);
		DVect2 JKR_Analytical(double overlap);

		// Contact model inheritance fields.
		quint32 inheritanceField_;

		// Effective translational stiffness.
		DVect2  effectiveTranslationalStiffness_;
		DAVect  effectiveRotationalStiffness_;      // (Twisting,Bending,Bending) Rotational stiffness (twisting always 0)

		// JKR model properties 
		double      shear_;			// Shear modulus
		double      poiss_;			// Poisson's ratio
		double      fric_;			// Coulomb friction coefficient
		DVect       jkr_F_;			// Force carried in the eepa model
		bool        jkr_S_;			// The current slip state
		double      rgap_;			// Reference gap 
		dpProps *   dpProps_;		// The viscous properties
		double      kn_;			// normal stiffness
		double      ks_;			// shear stiffness
		double      ks_fac_;		// shear stiffness scaling factor
		double		surf_adh_;		// surface adhesion/cohesion energy (energy/area = J/m^2)	
		int			active_mode_;	// mode=0: no negtaive overlap is allowed - contact activated and deactivated when physical contact is made (Simplified SJKR-A model)
									// mode=1: (default) negative overlap is allowed - contact activated on physical contact, but negative overlap (up to tear-off distance) allowed when unloading (full JKR model)
		double		distance_active_;// active distance - initially set to zero to establish contact, then after first time active it is set to the tear_off_ distance to allow for a seperation distance while under adhesion force
		double		a0_;			// contact patch radius where the JKR force is equal to zero
		double		pull_off_f_;	// the pull-off force which is the maximum tensile (adhesive) force
		double		tear_off_d_;	// tear-off distance where the the adhesion snaps to zero under te
		double      r_hertz_;		// effective contact radius: r_hertz_ = R1xR2/(R1+R2)
		double      c1_;			// constant used in the analytical solution to the JKR contact patch radius
												
		// rolling resistance properties
		double res_fric_;			// rolling friction coefficient
		DAVect res_M_;				// moment (bending only)         
		bool   res_S_;				// The current rolling resistance slip state
		double kr_;					// bending rotational stiffness (read-only, calculated internaly) 
		double fr_;					// rolling friction coefficient (rbar*res_fric_) (calculated internaly, not a property)
		double rbar_square_;		// used to update the rolling stiffness

		Energies *   energies_;		// The energies
	};
}
// EoF

Top

contactmodeljkr.cpp

// contactmodeljkr.cpp 
#include "contactmodeljkr.h"

#include "module/interface/icontactmechanical.h"
#include "module/interface/icontact.h"
#include "module/interface/ipiecemechanical.h"
#include "module/interface/ipiece.h"
#include "module/interface/ifishcalllist.h"

#include "../version.txt"

#include "utility/src/tptr.h"
#include "shared/src/mathutil.h"

#include "kernel/interface/iprogram.h"
#include "module/interface/icontactthermal.h"
#include "contactmodel/src/contactmodelthermal.h"

#include "module/interface/itablelist.h"
#include "module/interface/itable.h"

#include "utility/interface/ikerneloutput.h"
#include "utility/interface/itextoutput.h"
#include "fish/src/parameter.h"

#ifdef JKR_LIB
int __stdcall DllMain(void *, unsigned, void *) {
	return 1;
}

extern "C" EXPORT_TAG const char *getName() {
#if DIM==3
	return "contactmodelmechanical3dJKR";
#else
	return "contactmodelmechanical2dJKR";
#endif
}

extern "C" EXPORT_TAG unsigned getMajorVersion() {
	return MAJOR_VERSION;
}

extern "C" EXPORT_TAG unsigned getMinorVersion() {
	return MINOR_VERSION;
}

extern "C" EXPORT_TAG void *createInstance() {
	cmodelsxd::ContactModelJKR *m = new cmodelsxd::ContactModelJKR();
	return (void *)m;
}
#endif 

namespace cmodelsxd {
	static const quint32 shearMask = 0x00000002; // Base 1!
	static const quint32 poissMask = 0x00000004;
	static const quint32 fricMask = 0x00000008;
	static const quint32 resFricMask = 0x00004000;
	static const quint32 surfAdhMask = 0x00004000;

	using namespace itasca;

	int ContactModelJKR::index_ = -1;
	UInt ContactModelJKR::getMinorVersion() const { return MINOR_VERSION; }

	ContactModelJKR::ContactModelJKR() : inheritanceField_(shearMask | poissMask | fricMask | resFricMask | surfAdhMask)
		, effectiveTranslationalStiffness_(DVect2(0.0))
		, effectiveRotationalStiffness_(DAVect(0.0))
		, shear_(0.0)
		, poiss_(0.0)
		, fric_(0.0)
		, jkr_F_(DVect(0.0))
		, jkr_S_(false)
		, rgap_(0.0)
		, dpProps_(0)
		, res_fric_(0.0)
		, res_M_(DAVect(0.0))
		, res_S_(false)
		, kr_(0.0)
		, fr_(0.0)
		, kn_(0.0)			
		, ks_(0.0)					
		, ks_fac_(1.0)
		, surf_adh_(0.0)
		, active_mode_(1)
		, distance_active_(0.0)
		, a0_(0.0)
		, pull_off_f_(0.0)
		, tear_off_d_(0.0)
		, r_hertz_(0.0)
		, c1_(0.0)
		, rbar_square_(0.0)
		, energies_(0) {
	}

	ContactModelJKR::~ContactModelJKR() {
		// Make sure to clean up after yourself!
		if (dpProps_)
			delete dpProps_;
		if (energies_)
			delete energies_;
	}

	void ContactModelJKR::archive(ArchiveStream &stream) {
		// The stream allows one to archive the values of the contact model
		// so that it can be saved and restored. The minor version can be
		// used here to allow for incremental changes to the contact model too. 
		stream & shear_;
		stream & poiss_;
		stream & fric_;
		stream & jkr_F_;
		stream & jkr_S_;
		stream & rgap_;
		stream & res_fric_;
		stream & res_M_;
		stream & res_S_;
		stream & kr_;
		stream & fr_;
		stream & rbar_square_;
		stream & kn_;			
		stream & ks_;										
		stream & ks_fac_;
		stream & surf_adh_;
		stream & active_mode_;
		stream & distance_active_;
		stream & a0_;
		stream & pull_off_f_;
		stream & tear_off_d_;
		stream & r_hertz_;
		stream & c1_;
				
		if (stream.getArchiveState() == ArchiveStream::Save) {
			bool b = false;
			if (dpProps_) {
				b = true;
				stream & b;
				stream & dpProps_->dp_nratio_;
				stream & dpProps_->dp_sratio_;
				stream & dpProps_->dp_mode_;
				stream & dpProps_->dp_F_;
			}
			else
				stream & b;

			b = false;
			if (energies_) {
				b = true;
				stream & b;
				stream & energies_->estrain_;
				stream & energies_->errstrain_;
				stream & energies_->eslip_;
				stream & energies_->errslip_;
				stream & energies_->edashpot_;
			}
			else
				stream & b;
		}
		else {
			bool b(false);
			stream & b;
			if (b) {
				if (!dpProps_)
					dpProps_ = NEWC(dpProps());
				stream & dpProps_->dp_nratio_;
				stream & dpProps_->dp_sratio_;
				stream & dpProps_->dp_mode_;
				stream & dpProps_->dp_F_;
			}
			stream & b;
			if (b) {
				if (!energies_)
					energies_ = NEWC(Energies());
				stream & energies_->estrain_;
				stream & energies_->errstrain_;
				stream & energies_->eslip_;
				stream & energies_->errslip_;
				stream & energies_->edashpot_;
			}
		}

		stream & inheritanceField_;
		stream & effectiveTranslationalStiffness_;
		stream & effectiveRotationalStiffness_;
	}

	void ContactModelJKR::copy(const ContactModel *cm) {
		// Copy all of the contact model properties. Used in the CMAT 
		// when a new contact is created. 
		ContactModelMechanical::copy(cm);
		const ContactModelJKR *in = dynamic_cast<const ContactModelJKR*>(cm);
		if (!in) throw std::runtime_error("Internal error: contact model dynamic cast failed.");
		shear(in->shear());
		poiss(in->poiss());
		fric(in->fric());
		jkr_F(in->jkr_F());
		jkr_S(in->jkr_S());
		rgap(in->rgap());
		res_fric(in->res_fric());
		res_M(in->res_M());
		res_S(in->res_S());
		kr(in->kr());
		fr(in->fr());
		rbar_square(in->rbar_square());
		kn(in->kn());			
		ks(in->ks());	
		ks_fac(in->ks_fac());
		surf_adh(in->surf_adh());
		active_mode(in->active_mode());
		distance_active(in->distance_active());
		a0(in->a0());
		pull_off(in->pull_off());
		tear_off(in->tear_off());
		r_hertz(in->r_hertz());
		c1(in->c1());
				
		if (in->hasDamping()) {
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dp_nratio(in->dp_nratio());
			dp_sratio(in->dp_sratio());
			dp_mode(in->dp_mode());
			dp_F(in->dp_F());
		}
		if (in->hasEnergies()) {
			if (!energies_)
				energies_ = NEWC(Energies());
			estrain(in->estrain());
			errstrain(in->errstrain());
			eslip(in->eslip());
			errslip(in->errslip());
			edashpot(in->edashpot());
		}
		inheritanceField(in->inheritanceField());
		effectiveTranslationalStiffness(in->effectiveTranslationalStiffness());
		effectiveRotationalStiffness(in->effectiveRotationalStiffness());
	}

	QVariant ContactModelJKR::getProperty(uint i, const IContact *) const {
		// Return the property. The IContact pointer is provided so that 
		// more complicated properties, depending on contact characteristics,
		// can be calcualted. 
		QVariant var;
		switch (i) {
		case kwShear:     return shear_;
		case kwPoiss:     return poiss_;
		case kwFric:      return fric_;
		case kwjkrF:   var.setValue(jkr_F_); return var;
		case kwjkrS:   return jkr_S_;
		case kwDpNRatio:  return dpProps_ ? dpProps_->dp_nratio_ : 0;
		case kwDpSRatio:  return dpProps_ ? dpProps_->dp_sratio_ : 0;
		case kwDpMode:    return dpProps_ ? dpProps_->dp_mode_ : 0;
		case kwDpF: {
			dpProps_ ? var.setValue(dpProps_->dp_F_) : var.setValue(DVect(0.0));
			return var;
		}
		case kwResFric:    return res_fric_;
		case kwResMoment:  var.setValue(res_M_); return var;
		case kwResS:       return res_S_;
		
		case kwKsFac:		return ks_fac_;
		case kwSurfAdh:		return surf_adh_;
		case kwActiveMode:	return active_mode_;
		case kwA0:			return a0_;
		case kwPullOff:		return pull_off_f_;
		case kwTearOff:		return tear_off_d_;
		}
		assert(0);
		return QVariant();
	}

	bool ContactModelJKR::getPropertyGlobal(uint i) const {
		// Returns whether or not a property is held in the global axis system (TRUE)
		// or the local system (FALSE). Used by the plotting logic.
		switch (i) {
		case kwjkrF:
		case kwDpF:
		case kwResMoment:
			return false;
		}
		return true;
	}

	bool ContactModelJKR::setProperty(uint i, const QVariant &v, IContact *) {
		// Set a contact model property. Return value indicates that the timestep
		// should be recalculated. 
		dpProps dp;
		switch (i) {
		case kwShear: {
			if (!v.canConvert<double>())
				throw Exception("shear must be a double.");
			double val(v.toDouble());
			if (val <= 0.0)
				throw Exception("zero or negative shear not allowed in the 'full_jkr' model.");
			shear_ = val;
			return true;
		}
		case kwPoiss: {
			if (!v.canConvert<double>())
				throw Exception("poiss must be a double.");
			double val(v.toDouble());
			if (val < 0.0 || val > 0.5)
				throw Exception("poiss must be in the range [0, 0.5].");
			poiss_ = val;
			return true;
		}
		case kwFric: {
			if (!v.canConvert<double>())
				throw Exception("fric must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative fric not allowed.");
			fric_ = val;
			return false;
		}
		case kwjkrF: {
			if (!v.canConvert<DVect>())
				throw Exception("jkr_force must be a vector.");
			DVect val(v.value<DVect>());
			jkr_F_ = val;
			return false;
		}
		case kwDpNRatio: {
			if (!v.canConvert<double>())
				throw Exception("dp_nratio must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative dp_nratio not allowed.");
			if (val == 0.0 && !dpProps_)
				return false;
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_nratio_ = val;
			return true;
		}
		case kwDpSRatio: {
			if (!v.canConvert<double>())
				throw Exception("dp_sratio must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative dp_sratio not allowed.");
			if (val == 0.0 && !dpProps_)
				return false;
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_sratio_ = val;
			return true;
		}
		case kwDpMode: {
			if (!v.canConvert<int>())
				throw Exception("the viscous mode dp_mode must be 0 or 1");
			int val(v.toInt());
			if (val == 0 && !dpProps_)
				return false;
			if (val < 0 || val > 1)
				throw Exception("the viscous mode dp_mode must be 0 or 1.");
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_mode_ = val;
			return false;
		}
		case kwDpF: {
			if (!v.canConvert<DVect>())
				throw Exception("dp_force must be a vector.");
			DVect val(v.value<DVect>());
			if (val.fsum() == 0.0 && !dpProps_)
				return false;
			if (!dpProps_)
				dpProps_ = NEWC(dpProps());
			dpProps_->dp_F_ = val;
			return false;
		}
		case kwResFric: {
			if (!v.canConvert<double>())
				throw Exception("res_fric must be a double.");
			double val(v.toDouble());
			if (val < 0.0)
				throw Exception("negative res_fric not allowed.");
			res_fric_ = val;
			return true;
		}
		case kwResMoment: {
			DAVect val(0.0);
#ifdef TWOD               
			if (!v.canConvert<DAVect>() && !v.canConvert<double>())
				throw Exception("res_moment must be an angular vector.");
			if (v.canConvert<DAVect>())
				val = DAVect(v.value<DAVect>());
			else
				val = DAVect(v.toDouble());
#else
			if (!v.canConvert<DAVect>() && !v.canConvert<DVect>())
				throw Exception("res_moment must be an angular vector.");
			if (v.canConvert<DAVect>())
				val = DAVect(v.value<DAVect>());
			else
				val = DAVect(v.value<DVect>());
#endif
			res_M_ = val;
			return false;
		}		
		case kwKsFac: {
			if (!v.canConvert<double>())
				throw Exception("shear stiffness factor must be a double.");
			double val(v.toDouble());
			if (val <= 0.0)
				throw Exception("negative or zero ks_fac not allowed.");
			ks_fac_ = val;
			return true;
		}																								
		case kwSurfAdh: {
			if (!v.canConvert<double>())
				throw Exception("surface adhesion must be a double.");
			double val(v.toDouble());
			if (val <= 0.0)
				throw Exception("negative or zero surf_adh not allowed.");
			surf_adh_ = val;
			return true;
		}
		case kwActiveMode: {
			if (!v.canConvert<int>())
				throw Exception("the active mode must be an integer: 0 or 1.");
			int val(v.toInt());
			if (val < 0 || val > 1)
				throw Exception("active_mode must be 0 or 1.");
			active_mode_ = val;
			return true;
		}
		}//switch
		return false;
	}

	bool ContactModelJKR::getPropertyReadOnly(uint i) const {
		// Returns TRUE if a property is read only or FALSE otherwise. 
		switch (i) {
		case kwDpF:
		case kwjkrS:
		case kwResS:
		case kwA0:
		case kwPullOff:
		case kwTearOff:
			return true;
		default:
			break;
		}
		return false;
	}

	bool ContactModelJKR::supportsInheritance(uint i) const {
		// Returns TRUE if a property supports inheritance or FALSE otherwise. 
		switch (i) {
		case kwShear:
		case kwPoiss:
		case kwFric:
		case kwResFric:
		case kwSurfAdh:
			return true;
		default:
			break;
		}
		return false;
	}

	double ContactModelJKR::getEnergy(uint i) const {
		// Return an energy value. 
		double ret(0.0);
		if (!energies_)
			return ret;
		switch (i) {
		case kwEStrain:    return energies_->estrain_;
		case kwERRStrain:  return energies_->errstrain_;
		case kwESlip:      return energies_->eslip_;
		case kwERRSlip:    return energies_->errslip_;
		case kwEDashpot:   return energies_->edashpot_;
		}
		assert(0);
		return ret;
	}

	bool ContactModelJKR::getEnergyAccumulate(uint i) const {
		// Returns TRUE if the corresponding energy is accumulated or FALSE otherwise.
		switch (i) {
		case kwEStrain:   return true;
		case kwERRStrain: return false;
		case kwESlip:     return true;
		case kwERRSlip:   return true;
		case kwEDashpot:  return true;
		}
		assert(0);
		return false;
	}

	void ContactModelJKR::setEnergy(uint i, const double &d) {
		// Set an energy value. 
		if (!energies_) return;
		switch (i) {
		case kwEStrain:    energies_->estrain_ = d;		return;
		case kwERRStrain:  energies_->errstrain_ = d;	return;
		case kwESlip:      energies_->eslip_ = d;		return;
		case kwERRSlip:    energies_->errslip_ = d;		return;
		case kwEDashpot:   energies_->edashpot_ = d;	return;
		}
		assert(0);
		return;
	}

	bool ContactModelJKR::validate(ContactModelMechanicalState *state, const double &) {
		// Validate the / Prepare for entry into ForceDispLaw. The validate function is called when:
		// (1) the contact is created, (2) a property of the contact that returns a true via
		// the setProperty method has been modified and (3) when a set of cycles is executed
		// via the {cycle N} command.
		// Return value indicates contact activity (TRUE: active, FALSE: inactive).
		assert(state);
		const IContactMechanical *c = state->getMechanicalContact();
		assert(c);
		//
		if (state->trackEnergy_)
			activateEnergy();
		//		
		if ((inheritanceField_ & shearMask) || (inheritanceField_ & poissMask))
			updateEndStiffCoef(c);
		if (inheritanceField_ & fricMask)
			updateEndFric(c);
		if (inheritanceField_ & resFricMask)
			updateEndResFric(c);
		if (inheritanceField_ & surfAdhMask)
			updateEndSurfAdh(c);
		//
		if (shear_ <= 0.0) {
			throw Exception("'shear' must be specified using a value larger than zero in the 'full_jkr' model.");
		}
		if (surf_adh_ <= 0.0) {
			throw Exception("'surf_adh' must be specified using a value larger than zero in the 'full_jkr' model. For zero surface adhesion, use the Hertz model!");
		}
		//
		updateStiffCoef(c);								// calculate the stiffness values based on the material properties - specified directly or via inheritance
		updateEffectiveStiffness(state);				// effective stiffness for translation and rotation used in the time step estimation
		//
		return checkActivity(state->gap_);
	}

	void ContactModelJKR::updateStiffCoef(const IContactMechanical *con) {
		//
		double c12 = con->getEnd1Curvature().y();
		double c22 = con->getEnd2Curvature().y();
		r_hertz_ = c12 + c22;													
		if (r_hertz_ == 0.0)
			throw Exception("jkr contact model undefined for 2 non-curved surfaces");
		//
		r_hertz_ = 1.0 / r_hertz_;												// r_hertz = R1*R2/(R1 + R2)
		double young_eff = shear_ / (1.0 - poiss_);								// effective contact Young's modulus assuming the two contact pieces are identical in properties
		double shear_eff = shear_ / (4.0 - 2.0*poiss_);							// effective shear modulus for contact assuming the two contact pieces are identical in properties
		kn_ = 2.0 * young_eff;													// normal stiffness - needs to be multiplied by contact patch radius to obtain the tangent normal stiffness - (Marshall, 2009)
		ks_ = ks_fac_ * 8.0*shear_eff;											// shear stiffness - needs to be multiplied by contact patch radius to obtain the tangent shear stiffness - (Marshall, 2009)
		//
		a0_ = 9.0*M_PI*surf_adh_*r_hertz_*r_hertz_ / young_eff;					// contact patch radius where the JKR force is equal to zero
		a0_ = pow(a0_, 1.0 / 3.0);
		//
		pull_off_f_ = 3.0*M_PI*surf_adh_*r_hertz_;								// pull-off force, the maximum tensile/adhesion force in JKR model (not the force at tear-off when the contact snaps)
		tear_off_d_ = a0_* a0_ / (2.0*pow(6.0, (1.0 / 3.0))*r_hertz_);			// tear-off distance where the the adhesion snaps to zero under tension
		//
		c1_ = -4.0*M_PI*surf_adh_ * r_hertz_*r_hertz_ / young_eff;				// constant used in the analytical solution for the patch radius
		//
		// rolling resistance 
		kr_ = 0.0;
		fr_ = 0.0;
		if (res_fric_ > 0.0) {
			rbar_square_ = r_hertz_ * r_hertz_;			// store this value with contact - used again when 'ks_tangent' is updated during loading-unloading to calculate the rolling stiffness
			fr_ = res_fric_ * r_hertz_;					// store this value with contact - used in force-displacement calculation
			kr_ = ks_ * rbar_square_;					// based on 'ks_' for now, this is updated and based on the shear tangent stiffness during loading-unloading
		}
	}

	bool ContactModelJKR::endPropertyUpdated(const QString &name, const IContactMechanical *c) {
		// The endPropertyUpdated method is called whenever a surface property (with a name
		// that matches an inheritable contact model property name) of one of the contacting
		// pieces is modified. This allows the contact model to update its associated
		// properties. The return value denotes whether or not the update has affected
		// the time step computation (by having modified the translational or rotational
		// tangent stiffnesses). If true is returned, then the time step will be recomputed.  
		assert(c);
		QStringList availableProperties = getProperties().simplified().replace(" ", "").split(",", QString::SkipEmptyParts);
		QRegExp rx(name, Qt::CaseInsensitive);
		int idx = availableProperties.indexOf(rx) + 1;
		bool ret = false;

		if (idx <= 0)
			return ret;

		switch (idx) {
		case kwShear: { //Shear
			if (inheritanceField_ & shearMask)
				ret = true;							// return 'true' to ensure 'validate()' is called to affect the change in Shear Modulus through 'updateEndStiffCoef()'
			break;
		}
		case kwPoiss: { //Poisson's
			if (inheritanceField_ & poissMask)
				ret =true;							// return 'true' to ensure 'validate()' is called to affect the change in Poisson's ratio through 'updateEndStiffCoef()'
			break;
		}
		case kwFric: { //fric
			if (inheritanceField_ & fricMask)
				updateEndFric(c);					// friction does not influence any of the other parameters or time step size, so directly update the contact model friction
			break;
		}
		case kwResFric: { //rr_fric
			if (inheritanceField_ & resFricMask)
				ret = true;							// return 'true' to ensure 'validate()' is called to affect the change in rolling friction through 'updateEndResFric()' and 'fr_'
			break;
		}
		case kwSurfAdh: { //surf_adh
			if (inheritanceField_ & surfAdhMask)
				ret = true;							// return 'true' to ensure 'validate()' is called to affect the change in surface adhesion through 'updateEndSurfAdh()'
			break;
		}
	    //
		}
		return ret;
	}

	static const QString gstr("jkr_shear");
	static const QString nustr("jkr_poiss");
	bool ContactModelJKR::updateEndStiffCoef(const IContactMechanical *con) {
		assert(con);
		double g1 = shear_;
		double g2 = shear_;
		double nu1 = poiss_;
		double nu2 = poiss_;
		QVariant vg1 = con->getEnd1()->getProperty(gstr);
		QVariant vg2 = con->getEnd2()->getProperty(gstr);
		QVariant vnu1 = con->getEnd1()->getProperty(nustr);
		QVariant vnu2 = con->getEnd2()->getProperty(nustr);
		if (vg1.isValid() && vg2.isValid()) {
			g1 = vg1.toDouble();
			g2 = vg2.toDouble();
			if (g1 <= 0.0 || g2 <= 0.0)
				throw Exception("Negative or zero shear modulus not allowed in jkr contact model");
		}
		if (vnu1.isValid() && vnu2.isValid()) {
			nu1 = vnu1.toDouble();
			nu2 = vnu2.toDouble();
			if (nu1 < 0.0 || nu1 > 0.5 || nu2 < 0.0 || nu2 > 0.5)
				throw Exception("Poisson ratio should be in range [0,0.5] in jkr contact model");
		}
		if (g1*g2 == 0.0) return false;
		double es = 1.0 / ((1.0 - nu1) / (2.0*g1) + (1.0 - nu2) / (2.0*g2));
		double gs = 1.0 / ((2.0 - nu1) / g1 + (2.0 - nu2) / g2);
		poiss_ = (4.0*gs - es) / (2.0*gs - es);
		shear_ = 2.0*gs*(2 - poiss_);
		if (shear_ <= 0.0)
			throw Exception("Negative or zero shear modulus not allowed in jkr contact model");
		if (poiss_ < 0.0 || poiss_ > 0.5)
			throw Exception("Poisson ratio should be in range [0,0.5] in jkr contact model");
		return true;
	}
	
	static const QString fricstr("fric");
	bool ContactModelJKR::updateEndFric(const IContactMechanical *con) {
		assert(con);
		QVariant v1 = con->getEnd1()->getProperty(fricstr);
		QVariant v2 = con->getEnd2()->getProperty(fricstr);
		if (!v1.isValid() || !v2.isValid())
			return false;
		double fric1 = std::max(0.0, v1.toDouble());
		double fric2 = std::max(0.0, v2.toDouble());
		double val = fric_;
		fric_ = std::min(fric1, fric2);
		if (fric_ < 0.0)
			throw Exception("Negative friction value not allowed in jkr contact model");
		return ((fric_ != val));
	}

	static const QString rfricstr("rr_fric");
	bool ContactModelJKR::updateEndResFric(const IContactMechanical *con) {
		assert(con);
		QVariant v1 = con->getEnd1()->getProperty(rfricstr);
		QVariant v2 = con->getEnd2()->getProperty(rfricstr);
		if (!v1.isValid() || !v2.isValid())
			return false;
		double rfric1 = std::max(0.0, v1.toDouble());
		double rfric2 = std::max(0.0, v2.toDouble());
		double val = res_fric_;
		res_fric_ = std::min(rfric1, rfric2);
		if (res_fric_ < 0.0)
			throw Exception("Negative rolling friction value not allowed in jkr contact model");
		return ((res_fric_ != val));
	}

	static const QString surfadhstr("surf_adh");
	bool ContactModelJKR::updateEndSurfAdh(const IContactMechanical *con) {
		assert(con);
		QVariant v1 = con->getEnd1()->getProperty(surfadhstr);
		QVariant v2 = con->getEnd2()->getProperty(surfadhstr);
		if (!v1.isValid() || !v2.isValid())
			return false;
		double surfadh1 = std::max(0.0, v1.toDouble());
		double surfadh2 = std::max(0.0, v2.toDouble());
		if (surfadh1 <= 0.0 || surfadh2 <= 0.0)
			throw Exception("Negative or zero surface adhesion not allowed in 'full_jkr' contact model");
		double val = surf_adh_;
		surf_adh_ = 0.5*(surfadh1 + surfadh2);		// new surface adhesion enery
		return ((surf_adh_ != val));
	}

	void ContactModelJKR::updateEffectiveStiffness(ContactModelMechanicalState *state) {
		
		double overlap = rgap_ - state->gap_;
		if (overlap <= -distance_active_) {									// if contact was just created but not yet activated, distance_active_ == 0. 
			// The assumption below is based on the Hill contact model in PFC
			// 0.01% of Diameter = 0.0001*D			
			// For monodisperse system, r_hertz_ = 0.5*R where R is the particle radius
			// r_hertz_ = 0.5*(0.5*D) = 0.25*D where D is the particle diameter
			// D = 4*r_hertz_
			// 0.01%*D = 0.0001*D = 0.0001*4*r_hertz_ = 0.0004*r_hertz_
			overlap = 0.0004*r_hertz_;        
		}						
		//
		DVect2 force_radius = JKR_Analytical(overlap);						// analytical patch radius and force for this (assumed) overlap
		double patch_radius = force_radius.dof(1);							// contact patch radius
		// Assume Hertz tangent stiffness in the normal direction - the tangent stiffness for JKR is different, but no closed-form exists - good enough of an assumption
		double kn_tangent = kn_ * patch_radius;								// tangent stiffness in normal direction
		double ks_tangent = ks_ * patch_radius;								// tangent stiffness in the shear direction - tis is not an assumption, but the correct solution for JKR
		effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent);	// set using the tangent stiffness
		if (res_fric_ > 0.0) {												// rolling
			kr_ = ks_tangent * rbar_square_;								// based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
			effectiveRotationalStiffness_ = DAVect(kr_);					// effective rotational stiffness (bending only)
#if DIM==3 
			effectiveRotationalStiffness_.rx() = 0.0;
#endif
		}
		// correction if viscous damping active
		if (dpProps_) {
			DVect2 correct(1.0);
			if (dpProps_->dp_nratio_)
				correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
			if (dpProps_->dp_sratio_)
				correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
			effectiveTranslationalStiffness_ /= (correct*correct);
		}
	}
	//
	// --------------------------------------------------------------------------------------------------------------------------------------------------------
    //  
	DVect2 ContactModelJKR::JKR_Analytical(double overlap) {															// analytical solution according to Parteli et al. 2014
		double rad_over = r_hertz_ * overlap;
		double c0 = rad_over * rad_over;											
		double c2 = -2.0*rad_over;													
		double P = -1.0 / 12.0*c2*c2 - c0;											
		double Q = -1.0 / 108.0*c2*c2*c2 + c0 * c2 / 3.0 - c1_*c1_ / 8.0;			
		double U = pow((-0.5*Q + sqrt(0.25*Q*Q + P*P*P / 27.0)), (1.0 / 3.0));		
		double s = 0.0;
		if (P == 0.0) {
			s = -5.0*c2 / 6.0 - pow(Q, (1.0 / 3.0));
		}
		else {
			s = -5.0*c2 / 6.0 + U - P / (3.0*U);
		}
		double w = sqrt(c2 + 2.0*s);
		double lambda = c1_ / (2.0*w);
		double patch_radius = 0.5*(w + sqrt(w*w - 4.0*(c2 + s + lambda)));			
		double patch_radius3 = patch_radius* patch_radius*patch_radius;				
		DVect2 ret(0.0);
		// force = 4*E*(a^3)/(3*R) - sqrt(16*pi*surf_adh*E*a^3);
		// with E = kn/2 - where E is the effective contact Young's modulus - see updateStiffCoef(const IContactMechanical *con)
		// the above reduces to the following for the JKR force ...
		ret.rdof(0) = (2.0 / 3.0)*kn_*patch_radius3 / r_hertz_ - sqrt(8.0*M_PI*surf_adh_*kn_*patch_radius3);			// set the total force
		ret.rdof(1) = patch_radius;																						// set the patch contact radius
		return ret;
	}
	//
	// --------------------------------------------------------------------------------------------------------------------------------------------------------
	//  
	bool ContactModelJKR::forceDisplacementLaw(ContactModelMechanicalState *state, const double &timestep) {
		assert(state);
		//
		// Current overlap
		double overlap = rgap_ - state->gap_;							// state->gap is negative in sign when the particles physically overlap. If rgap == 0, then the overlap is positive when the two particles physically overlap.
		// Relative translational increment
		DVect trans = state->relativeTranslationalIncrement_;
		//
		if (active_mode_ == 1 && distance_active_ == 0.0) {				// only when in 'full jkr' mode = 1, for Simplified JKR-A (SJKR-A), negative overlap is not allowed
			distance_active_ = tear_off_d_;								// contact has become active - set the activity distance to rgap_ + distance_active = rgap_ + tear_off to allow for negative overlap under adhesion force
		}
		//
		// The contact was just activated from an inactive state
		// Trigger the FISH callback if one is hooked up to the contact_activated event.
		if (cmEvents_[fActivated] >= 0) {
			// An FArray of QVariant is returned and these will be passed
			// to the FISH function as an array of FISH symbols as the second
			// argument to the FISH callback function. 
			auto c = state->getContact();
			std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()) };
			IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
			fi->setCMFishCallArguments(c, arg, cmEvents_[fActivated]);

		}
		//
		// Angular dispacement increment.
		DAVect ang = state->relativeAngularIncrement_;
		//
		DVect jkr_F_old = jkr_F_;													// store for energy computations only
		//
		// NORMAL FORCE ==============================================================================================================	
		//
		DVect2 force_radius = JKR_Analytical(overlap);								// analytical force and patch radius
		jkr_F_.rdof(0) = force_radius.dof(0);										// current contact force in the normal direction
		double patch_radius = force_radius.dof(1);									// current contact patch radius
		//
		double kn_tangent = std::abs((jkr_F_old.dof(0) - jkr_F_.dof(0)) / trans.dof(0));// tangent stiffness in normal direction - used for timestep calculations and damping
		//
		// SHEAR FORCE ==============================================================================================================
		//
		double ks_tangent = ks_* patch_radius;										// tangent stiffness in the shear direction
		DVect sforce(0.0);
		// dim holds the dimension (e.g., 2 for 2D and 3 for 3D)
		// Loop over the shear components (note: the 0 component is the normal component) and calculate the shear force.
		for (int i = 1; i < dim; ++i)
			sforce.rdof(i) = jkr_F_.dof(i) - trans.dof(i) * ks_tangent;				// shear force update
		//
		// The canFail flag corresponds to whether or not the contact can undergo non-linear
		// force-displacement response. If the SOLVE ELASTIC command is given then the 
		// canFail state is set to FALSE. Otherwise it is always TRUE. 
		if (state->canFail_) {
			// Resolve sliding. This is the normal force multiplied by the coefficient of friction.
			double crit = fric_ * std::abs(jkr_F_.rdof(0) + 2.0*pull_off_f_);								
			// The is the magnitude of the shear force.
			double sfmag = sforce.mag();
			// Sliding occurs when the magnitude of the shear force is greater than the critical value.
			if (sfmag > crit) {
				// Lower the shear force to the critical value for sliding.
				double rat = crit / sfmag;
				sforce *= rat;
				// Handle the slip_change event if one has been hooked up. Sliding has commenced.  
				if (!jkr_S_ && cmEvents_[fSlipChange] >= 0) {
					auto c = state->getContact();
					std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
														 fish::Parameter((qint64)0) };
					IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
					fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
				}
				jkr_S_ = true;
			}
			else {
				// Handle the slip_change event if one has been hooked up and
				// the contact was previously sliding. Sliding has ceased.  
				if (jkr_S_) {
					if (cmEvents_[fSlipChange] >= 0) {
						auto c = state->getContact();
						std::vector<fish::Parameter> arg = { fish::Parameter(c->getIThing()),
															 fish::Parameter((qint64)1) };
						IFishCallList* fi = const_cast<IFishCallList*>(state->getProgram()->findInterface<IFishCallList>());
						fi->setCMFishCallArguments(c, arg, cmEvents_[fSlipChange]);
					}
					jkr_S_ = false;
				}
			}
		}
		//
		// Set the shear components of the total force.
		for (int i = 1; i < dim; ++i)
			jkr_F_.rdof(i) = sforce.dof(i);
		//
		// ROLLING RESISTANCE ==============================================================================================================
		//
		DAVect res_M_old = res_M_;													// used in energy calculation only
		if (fr_ == 0.0) {											
			res_M_.fill(0.0);
			kr_ = 0.0;
		}
		else {
			DAVect angStiff(0.0);
			DAVect MomentInc(0.0);
			kr_ = ks_tangent * rbar_square_;										// update rolling stiffness based on the tangent shear stiffness (Wensrich and Katterfeld, 2012)
#if DIM==3 
			angStiff.rx() = 0.0;
			angStiff.ry() = kr_;
#endif
			angStiff.rz() = kr_;
			MomentInc = ang * angStiff * (-1.0);
			res_M_ += MomentInc;
			if (state->canFail_) {
				// Account for bending strength
				double rmax = std::abs(fr_*(jkr_F_.rdof(0) + 2.0*pull_off_f_));		// Using the same normal force as in the shear limit
				double rmag = res_M_.mag();
				if (rmag > rmax) {
					double fac = rmax / rmag;
					res_M_ *= fac;
					res_S_ = true;
				}
				else {
					res_S_ = false;
				}
			}
		}
		//
		// STIFFNESS FOR TIME STEP UPDATES ============================================================================================
		//
		// set effective stiffness in the normal and shear directions for timestep calculation
		effectiveTranslationalStiffness_ = DVect2(kn_tangent, ks_tangent);
		// set the effective rotational stiffness (bending only)
		effectiveRotationalStiffness_ = DAVect(kr_);
#if DIM==3 
		effectiveRotationalStiffness_.rx() = 0.0;
#endif
		//
		// DAMPING FORCE ==============================================================================================================
		//
		// Account for dashpot forces if the dashpot structure has been defined. 
		if (dpProps_) {
			dpProps_->dp_F_.fill(0.0);
			double vcn(0.0), vcs(0.0);
			// Calculate the damping coefficients. 
			setDampCoefficients(state->inertialMass_, &kn_tangent, &ks_tangent, &vcn, &vcs);				// this is based on the current tangent stiffness
			// First damp the shear components
			for (int i = 1; i < dim; ++i)
				dpProps_->dp_F_.rdof(i) = trans.dof(i) * (-1.0* vcs) / timestep;
			// Damp the normal component
			dpProps_->dp_F_.rx() -= trans.x() * vcn / timestep;
			//
			if (jkr_S_ && dpProps_->dp_mode_ == 1) {
				// Limit in shear if sliding.
				double dfn = dpProps_->dp_F_.rx();
				dpProps_->dp_F_.fill(0.0);
				dpProps_->dp_F_.rx() = dfn;
			}
			// Correct effective translational stiffness
			DVect2 correct(1.0);
			if (dpProps_->dp_nratio_)
				correct.rx() = sqrt(1.0 + dpProps_->dp_nratio_*dpProps_->dp_nratio_) - dpProps_->dp_nratio_;
			if (dpProps_->dp_sratio_)
				correct.ry() = sqrt(1.0 + dpProps_->dp_sratio_*dpProps_->dp_sratio_) - dpProps_->dp_sratio_;
			effectiveTranslationalStiffness_ /= (correct*correct);
		}
		
		//Compute energies if energy tracking has been enabled. 
		if (state->trackEnergy_) {
			assert(energies_);
			// calculate the strain energy increment in the normal direction				 
			energies_->estrain_ -= 0.5*(jkr_F_old.dof(0) + jkr_F_.dof(0))*trans.dof(0);// under loading, trans.x < 0, hence the negative sign
			if (ks_) {
				DVect u_s_elastic = trans;										// set the elastic displacement increment equal to the total displacement increment
				u_s_elastic.rx() = 0.0;											// set the normal component to zero: u_s_elatic = [0, trans_shear_1, trans_shear_2]
				DVect shearF = jkr_F_;											// set the shear force equal to the total jkr force (including the normal component)
				shearF.rx() = 0.0;												// set the normal component to zero: shearF = [0, shear_force_1, shear_force_2]
				jkr_F_old.rx() = 0.0;											// set normal component of the previous force equal to zero
				DVect avg_F_s = (shearF + jkr_F_old)*0.5;						// average shear force vector
				if (jkr_S_) {													// if sliding, calculate the slip energy and accumulate it
					DVect u_s_total = u_s_elastic;								// total shear displacement increment
					u_s_elastic = (shearF - jkr_F_old) / ks_tangent;			// elastic shear displacement increment
					energies_->eslip_ -= std::min(0.0, (avg_F_s | (u_s_total + u_s_elastic))); // where (u_s_total + u_s_elatic) is the slip displacment increment (due to the sign convention, the terms are added up)
				}
				energies_->estrain_ -= avg_F_s | u_s_elastic;					// add the elastic component (if any - if previous force equals current force, the elastic incrment is zero)
			}
			// Add the rolling resistance energy contributions.							// done incrementally since the stiffness kr_ changes with the tangent shear stiffness (non-linearly)
			if (kr_) {
				DAVect t_s_elastic = ang;												// set the elastic rotation increment equal to the total angle increment
				DAVect avg_M = (res_M_ + res_M_old)*0.5;								// average moment from this time step and the previous
				if (res_S_) {															// if sliding, calculate the slip energy and accumulate it
					t_s_elastic = (res_M_ - res_M_old) / kr_;							// elastic angle increment
					energies_->errslip_ -= std::min(0.0, (avg_M | (ang + t_s_elastic)));// where (ang + t_s_elatic) is the slip rotation increment (due to the sign convention, the terms are added up)
				}	
				energies_->errstrain_ -= avg_M | t_s_elastic;							// add the elastic component (if any - if previous force equals current force, the elastic incrment is zero)
			}
			// Calculate damping energy (accumulated) if the dashpots are active. 
			if (dpProps_) {
				energies_->edashpot_ -= dpProps_->dp_F_ | trans;
			}
		}

		// This is just a sanity check to ensure, in debug mode, that the force isn't wonky. 
		assert(jkr_F_ == jkr_F_);
		return true;
	}

	void ContactModelJKR::setForce(const DVect &, IContact *) {
		//
	}

	DVect ContactModelJKR::getForce(const IContactMechanical *) const {
		DVect ret(jkr_F_);
		if (dpProps_)
			ret += dpProps_->dp_F_;																		
		return ret;
	}

	DAVect ContactModelJKR::getMomentOn1(const IContactMechanical *c) const {
		DVect force = getForce(c);
		DAVect ret(res_M_);
		c->updateResultingTorqueOn1Local(force, &ret);
		return ret;
	}

	DAVect ContactModelJKR::getMomentOn2(const IContactMechanical *c) const {
		DVect force = getForce(c);
		DAVect ret(res_M_);
		c->updateResultingTorqueOn2Local(force, &ret);
		return ret;
	}

	void ContactModelJKR::setDampCoefficients(const double &mass, double *kn_tangent, double *ks_tangent, double *vcn, double *vcs) {
		*vcn = dpProps_->dp_nratio_ * 2.0 * sqrt(mass*(*kn_tangent));
		*vcs = dpProps_->dp_sratio_ * 2.0 * sqrt(mass*(*ks_tangent));
	}

} // namespace cmodelsxd
// EoF

Top