yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01167
[svn] r1730 - trunk/extra
Author: eudoxos
Date: 2009-03-25 22:30:06 +0100 (Wed, 25 Mar 2009)
New Revision: 1730
Modified:
trunk/extra/Brefcom.cpp
trunk/extra/Brefcom.hpp
Log:
1. Add rate-dependent damage to normal and viscoplasticity to shear components of Brefcom (not yet tested, just compiles); other cleanups there.
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-03-24 21:02:59 UTC (rev 1729)
+++ trunk/extra/Brefcom.cpp 2009-03-25 21:30:06 UTC (rev 1730)
@@ -110,8 +110,10 @@
if(neverDamage) contPhys->neverDamage=true;
if(cohesiveThresholdIter<0 || Omega::instance().getCurrentIteration()<cohesiveThresholdIter) contPhys->isCohesive=true;
else contPhys->isCohesive=false;
- contPhys->tau=tau;
- contPhys->expDmgRate=expDmgRate;
+ contPhys->dmgTau=dmgTau;
+ contPhys->dmgRateExp=plRateExp;
+ contPhys->plTau=plTau;
+ contPhys->plRateExp=plRateExp;
interaction->interactionPhysics=contPhys;
}
@@ -143,6 +145,40 @@
CREATE_LOGGER(ef2_Spheres_Brefcom_BrefcomLaw);
+Real BrefcomContact::solveBeta(const Real c, const Real N){
+ const int maxIter=20;
+ const Real maxError=1e-12;
+ Real f, ret=0.;
+ for(int i=0; i<maxIter; i++){
+ Real aux=c*exp(N*ret)+exp(ret);
+ f=log(aux);
+ if(fabs(f)<maxError) return ret;
+ Real df=(c*N*exp(N*ret)+exp(ret))/aux;
+ ret-=f/df;
+ }
+ LOG_FATAL("No convergence, c="<<c<<", ret="<<ret<<", f="<<f);
+ throw runtime_error("ef2_Spheres_Brefcom_BrefcomLaw::solveBeta failed to converge.");
+}
+
+Real BrefcomContact::computeDmgOverstress(Real epsN, Real dt){
+ if(kappaD>=epsN*omega){ // unloading, no viscous stress
+ kappaD=epsN*omega;
+ return 0.0;
+ }
+ Real c=epsCrackOnset*(1-omega)*pow(dmgTau/dt,dmgRateExp)*pow(epsN*omega-kappaD,dmgRateExp-1.);
+ Real beta=solveBeta(c,dmgRateExp);
+ Real deltaDmgStrain=(epsN*omega-kappaD)*exp(beta);
+ kappaD+=deltaDmgStrain;
+ return (epsN*omega-kappaD)*E;
+}
+
+Real BrefcomContact::computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt){
+ if(sigmaTNorm<sigmaTYield) return 1.;
+ Real c=/* should this be sigmaT0?? */ sigmaTNorm*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-1.);
+ Real beta=solveBeta(c,plRateExp);
+ return 1.-exp(beta)*(1-sigmaTYield/sigmaTNorm);
+}
+
void ef2_Spheres_Brefcom_BrefcomLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
//timingDeltas->start();
SpheresContactGeometry* contGeom=static_cast<SpheresContactGeometry*>(_geom.get());
@@ -155,7 +191,7 @@
if(BC->omega>=1.0 && BC->omegaThreshold>=1.0) return;
// shorthands
- Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& tau(BC->tau); const Real& expDmgRate(BC->expDmgRate); const Real& omegaThreshold(BC->omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& dt=Omega::instance().getTimeStep(); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage);
+ Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(BC->omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& dt=Omega::instance().getTimeStep(); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const Real& plTau(BC->plTau);
/* const Real& transStrainCoeff(BC->transStrainCoeff); const Real& epsTrans(BC->epsTrans); const Real& xiShear(BC->xiShear); */
Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
@@ -198,71 +234,11 @@
void BrefcomLaw::action(MetaBody* _rootBody){
rootBody=_rootBody;
- if(useFunctor){ // testing the functor
- if(!functor) functor=shared_ptr<ef2_Spheres_Brefcom_BrefcomLaw>(new ef2_Spheres_Brefcom_BrefcomLaw);
- functor->logStrain=logStrain;
- FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
- if(!I->isReal) continue;
- functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), rootBody);
- }
- return;
- }
-
+ if(!functor) functor=shared_ptr<ef2_Spheres_Brefcom_BrefcomLaw>(new ef2_Spheres_Brefcom_BrefcomLaw);
+ functor->logStrain=logStrain;
FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
if(!I->isReal) continue;
- BC=YADE_PTR_CAST<BrefcomContact>(I->interactionPhysics);
- contGeom=YADE_PTR_CAST<SpheresContactGeometry>(I->interactionGeometry);
- assert(BC); assert(contGeom);
- /* kept fully damaged contacts; note that normally the contact is deleted _after_ the BREFCOM_MATERIAL_MODEL,
- * i.e. if it is 1.0 here, omegaThreshold is >= 1.0 for sure.
- * &&'ing that just to make sure anyway ...
- */
- if(BC->omega>=1.0 && BC->omegaThreshold>=1.0) continue;
-
- // shorthands
- Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); /* const Real& xiShear(BC->xiShear);*/ const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& tau(BC->tau); const Real& expDmgRate(BC->expDmgRate); const Real& omegaThreshold(BC->omegaThreshold); /* const Real& transStrainCoeff(BC->transStrainCoeff); const Real& epsTrans(BC->epsTrans); */ const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage);
- // for python access
- Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs);
- // for rate-dependence
- const Real& dt=Omega::instance().getTimeStep();
-
- assert(contGeom->hasShear);
-
- epsN=contGeom->epsN();
- epsT=contGeom->epsT();
-
- contGeom->relocateContactPoints(); // allow very large mutual rotations
-
- if(logStrain && epsN<0){
- Real epsN0=epsN;
- epsN=log(epsN0+1);
- epsT*=epsN/epsN0;
- }
-
- #ifdef BREFCOM_MATERIAL_MODEL
- BREFCOM_MATERIAL_MODEL
- #else
- sigmaN=E*epsN;
- sigmaT=G*epsT;
- #endif
-
- if(omega>omegaThreshold){
- I->isReal=false;
- const shared_ptr<Body>& body1=Body::byId(I->getId1(),_rootBody), body2=Body::byId(I->getId2(),_rootBody); assert(body1); assert(body2);
- const shared_ptr<BrefcomPhysParams>& rbp1=YADE_PTR_CAST<BrefcomPhysParams>(body1->physicalParameters), rbp2=YADE_PTR_CAST<BrefcomPhysParams>(body2->physicalParameters);
- if(BC->isCohesive){rbp1->numBrokenCohesive+=1; rbp2->numBrokenCohesive+=1; rbp1->epsPlBroken+=epsPlSum; rbp2->epsPlBroken+=epsPlSum;}
- LOG_DEBUG("Contact #"<<I->getId1()<<"=#"<<I->getId2()<<" is damaged over thershold ("<<omega<<">"<<omegaThreshold<<") and has been deleted (isReal="<<I->isReal<<")");
- continue;
- }
-
- #define NNAN(a) assert(!isnan(a));
- #define NNANV(v) assert(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
- // store Fn (and Fs?), for use with GlobalStiffnessCounter?
- NNAN(sigmaN); NNANV(sigmaT);
- NNAN(crossSection);
- Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
- Fs=sigmaT*crossSection; BC->shearForce=Fs;
- applyForce(crossSection*(contGeom->normal*sigmaN + sigmaT),I->getId1(),I->getId2()); /* this is the force applied on the _first_ body; inverted applied to the second */
+ functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), rootBody);
}
}
Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp 2009-03-24 21:02:59 UTC (rev 1729)
+++ trunk/extra/Brefcom.hpp 2009-03-25 21:30:06 UTC (rev 1730)
@@ -70,10 +70,14 @@
omegaThreshold,
//! weight coefficient for shear strain when computing the strain semi-norm kappaD
xiShear,
- //! characteristic time (if non-positive, the law without rate-dependence is used)
- tau,
+ //! characteristic time for damage (if non-positive, the law without rate-dependence is used)
+ dmgTau,
//! exponent in the rate-dependent damage evolution
- expDmgRate,
+ dmgRateExp,
+ //! characteristic time for viscoplasticity (if non-positive, no rate-dependence for shear)
+ plTau,
+ //! exponent in the rate-dependent viscoplasticity
+ plRateExp,
//! coefficient that takes transversal strain into accound when calculating kappaDReduced
transStrainCoeff;
/*! Up to now maximum normal strain (semi-norm), non-decreasing in time. */
@@ -91,7 +95,14 @@
// FIXME: Fn and Fs are stored as Vector3r normalForce, shearForce in NormalShearInteraction
Real omega, Fn, sigmaN, epsN, relResidualStrength; Vector3r epsT, sigmaT, Fs;
- BrefcomContact(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), tau(0), expDmgRate(0), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; }
+
+ static Real solveBeta(const Real c, const Real N);
+ Real computeDmgOverstress(Real epsN,Real dt);
+ Real computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt);
+
+
+
+ BrefcomContact(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), dmgTau(-1), dmgRateExp(0), plTau(-1), plRateExp(0), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; }
// BrefcomContact(Real _E, Real _G, Real _tanFrictionAngle, Real _undamagedCohesion, Real _equilibriumDist, Real _crossSection, Real _epsCrackOnset, Real _epsFracture, Real _expBending, Real _xiShear, Real _tau=0, Real _expDmgRate=1): InteractionPhysics(), E(_E), G(_G), tanFrictionAngle(_tanFrictionAngle), undamagedCohesion(_undamagedCohesion), equilibriumDist(_equilibriumDist), crossSection(_crossSection), epsCrackOnset(_epsCrackOnset), epsFracture(_epsFracture), expBending(_expBending), xiShear(_xiShear), tau(_tau), expDmgRate(_expDmgRate) { epsT=Vector3r::ZERO; kappaD=0; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; /*TRVAR5(epsCrackOnset,epsFracture,Kn,crossSection,equilibriumDist); */ }
virtual ~BrefcomContact();
@@ -105,8 +116,10 @@
(epsFracture)
(omegaThreshold)
(xiShear)
- (tau)
- (expDmgRate)
+ (dmgTau)
+ (dmgRateExp)
+ (plTau)
+ (plRateExp)
(transStrainCoeff)
(kappaD)
@@ -153,14 +166,14 @@
REGISTER_SERIALIZABLE(BrefcomPhysParams);
class ef2_Spheres_Brefcom_BrefcomLaw: public ConstitutiveLaw{
+ public:
/*! Cohesion evolution law (it is 1-funcH here) */
Real funcH(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{ return 1-funcG(kappaD,epsCrackOnset,epsFracture,neverDamage); }
/*! Damage evolution law */
- inline Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{
+ static Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) {
if(kappaD<epsCrackOnset || neverDamage) return 0;
return 1.-(epsCrackOnset/kappaD)*exp(-(kappaD-epsCrackOnset)/epsFracture);
}
- public:
bool logStrain;
ef2_Spheres_Brefcom_BrefcomLaw(): logStrain(false){ /*timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/ }
void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
@@ -183,20 +196,16 @@
/*! Cohesion evolution law (it is 1-funcH here) */
Real funcH(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{ return 1-funcG(kappaD,epsCrackOnset,epsFracture,neverDamage); }
/*! Damage evolution law */
- inline Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{
- if(kappaD<epsCrackOnset || neverDamage) return 0;
- return 1.-(epsCrackOnset/kappaD)*exp(-(kappaD-epsCrackOnset)/epsFracture);
- }
+ Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) const{ return ef2_Spheres_Brefcom_BrefcomLaw::funcG(kappaD,epsCrackOnset,epsFracture,neverDamage); }
public:
bool logStrain;
- bool useFunctor;
- BrefcomLaw(): logStrain(false),useFunctor(false) { Shop::Bex::initCache(); };
+ BrefcomLaw(): logStrain(false) { Shop::Bex::initCache(); };
void action(MetaBody*);
protected:
NEEDS_BEX("Force","Momentum");
REGISTER_CLASS_AND_BASE(BrefcomLaw,InteractionSolver);
- REGISTER_ATTRIBUTES(InteractionSolver,(logStrain)(useFunctor));
+ REGISTER_ATTRIBUTES(InteractionSolver,(logStrain));
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(BrefcomLaw);
@@ -216,7 +225,7 @@
expBending is positive if the damage evolution function is concave after fracture onset;
reasonable value seems like 4.
*/
- Real sigmaT, expBending, xiShear, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, omegaThreshold, transStrainCoeff;
+ Real sigmaT, expBending, xiShear, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, omegaThreshold, transStrainCoeff, dmgTau, dmgRateExp, plTau, plRateExp;
//! Should new contacts be cohesive? They will before this iter#, they will not be afterwards. If 0, they will never be. If negative, they will always be created as cohesive.
long cohesiveThresholdIter;
//! Create contacts that don't receive any damage (BrefcomContact::neverDamage=true); defaults to false
@@ -228,8 +237,8 @@
xiShear=0;
neverDamage=false;
cohesiveThresholdIter=-1;
- tau=-1; expDmgRate=0;
- omegaThreshold=0.98;
+ dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
+ omegaThreshold=0.999;
}
virtual void go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction);
@@ -242,8 +251,10 @@
(neverDamage)
(epsCrackOnset)
(relDuctility)
- (tau)
- (expDmgRate)
+ (dmgTau)
+ (dmgRateExp)
+ (plTau)
+ (plRateExp)
(omegaThreshold)
(transStrainCoeff)
);
