yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3860
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Fri 2016-05-13 15:07:22 +0200
message:
  Modification of CPM
  
  added parameter to combine tensile and shear strain for damage evaluation
  added plasticity in compression
  fixed one example (one depricated parameter)
modified:
  examples/concrete/periodic.py
  pkg/dem/ConcretePM.cpp
  pkg/dem/ConcretePM.hpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'examples/concrete/periodic.py'
--- examples/concrete/periodic.py	2013-04-24 19:49:00 +0000
+++ examples/concrete/periodic.py	2016-05-13 13:07:22 +0000
@@ -39,7 +39,7 @@
 	sigmaT=3.5e6,
 	frictionAngle=atan(0.8),
 	epsCrackOnset=1e-4,
-	crackOpening=1e-6,
+	relDuctility = 30,
 
 	intRadius=1.5,
 	dtSafety=.4,
@@ -70,7 +70,7 @@
 # load the packing (again);
 #
 import cPickle as pickle
-concreteId=O.materials.append(CpmMat(young=young, frictionAngle=frictionAngle, density=4800, sigmaT=sigmaT, crackOpening=crackOpening, epsCrackOnset=epsCrackOnset, poisson=poisson, isoPrestress=isoPrestress))
+concreteId=O.materials.append(CpmMat(young=young, frictionAngle=frictionAngle, density=4800, sigmaT=sigmaT, relDuctility=relDuctility, epsCrackOnset=epsCrackOnset, poisson=poisson, isoPrestress=isoPrestress))
 sphDict=pickle.load(open(packingFile))
 from yade import pack
 sp=pack.SpherePack()

=== modified file 'pkg/dem/ConcretePM.cpp'
--- pkg/dem/ConcretePM.cpp	2016-04-23 17:54:37 +0000
+++ pkg/dem/ConcretePM.cpp	2016-05-13 13:07:22 +0000
@@ -66,6 +66,7 @@
 		#define _CPATTR(a) cpmPhys->a=mat1->a
 			_CPATTR(epsCrackOnset);
 			_CPATTR(relDuctility);
+			_CPATTR(equivStrainShearContrib);
 			_CPATTR(neverDamage);
 			_CPATTR(dmgTau);
 			_CPATTR(dmgRateExp);
@@ -84,6 +85,7 @@
 			cpmPhys->isCohesive = (cohesiveThresholdIter < 0 || scene->iter < cohesiveThresholdIter);
 			_AVGATTR(epsCrackOnset);
 			_AVGATTR(relDuctility);
+			_AVGATTR(equivStrainShearContrib);
 			cpmPhys->neverDamage = (mat1->neverDamage || mat2->neverDamage);
 			_AVGATTR(dmgTau);
 			_AVGATTR(dmgRateExp);
@@ -307,7 +309,7 @@
 		if (b1index == sphereIndex && b2index == sphereIndex) { // both bodies are spheres
 			const Vector3r& pos1 = Body::byId(I->id1,scene)->state->pos;
 			const Vector3r& pos2 = Body::byId(I->id2,scene)->state->pos;
-			Real minRad = (geom->refR1 <= 0? geom->refR2 : (geom->refR2 <=0? geom->refR1 : min(geom->refR1,geom->refR2)));
+			Real minRad = (geom->refR1 <= 0? geom->refR2 : (geom->refR2 <=0? geom->refR1 : std::min(geom->refR1,geom->refR2)));
 			Vector3r shift2 = scene->isPeriodic? Vector3r(scene->cell->hSize*I->cellDist.cast<Real>()) : Vector3r::Zero();
 			phys->refLength = (pos2 - pos1 + shift2).norm();
 			phys->crossSection = Mathr::PI*pow(minRad,2);
@@ -353,16 +355,16 @@
 
 	#ifdef CPM_MATERIAL_MODEL
 		Vector3r& epsTPl(phys->epsTPl);
-		Real& epsNPl(phys->epsNPl);
 		const Real& dt = scene->dt;
 		const Real& dmgTau(phys->dmgTau);
 		const Real& plTau(phys->plTau);
 		const Real& yieldLogSpeed(this->yieldLogSpeed);
 		const int& yieldSurfType(this->yieldSurfType);
 		const Real& yieldEllipseShift(this->yieldEllipseShift);
-		const Real& epsSoft(this->epsSoft);
-		const Real& relKnSoft(this->relKnSoft);
 	#endif
+	Real& epsNPl(phys->epsNPl);
+	const Real& epsSoft(this->epsSoft);
+	const Real& relKnSoft(this->relKnSoft);
 
 	TIMING_DELTAS_CHECKPOINT("GO A");
 	
@@ -384,11 +386,17 @@
 		/* simplified public model */
 		epsN += phys->isoPrestress/E;
 		/* very simplified version of the constitutive law */
-		kappaD = max(max((Real)0.,epsN),kappaD); /* internal variable, max positive strain (non-decreasing) */
+		Real xi2 = std::pow(phys->equivStrainShearContrib,2);
+		Real epsNorm = std::sqrt(std::pow(std::max(epsN-epsNPl,0.),2)+xi2*epsT.squaredNorm());
+		kappaD = std::max(epsNorm,kappaD); /* internal variable, max positive strain (non-decreasing) */
 		omega = isCohesive? phys->funcG(kappaD,epsCrackOnset,epsFracture,neverDamage,damLaw) : 1.; /* damage variable (non-decreasing, as funcG is also non-decreasing) */
-		sigmaN = (1-(epsN>0?omega:0))*E*epsN; /* damage taken in account in tension only */
+		sigmaN = (1-(epsN-epsNPl>0?omega:0))*E*(epsN-epsNPl); /* damage taken in account in tension only */
+		if((epsSoft<0) && (epsN-epsNPl<epsSoft)){ /* plastic slip in compression */
+			Real sigmaNSoft=E*(epsSoft+relKnSoft*(epsN-epsNPl-epsSoft));
+			if(sigmaNSoft>sigmaN){ /*assert(sigmaNSoft>sigmaN);*/ epsNPl+=(sigmaN-sigmaNSoft)/E; sigmaN=sigmaNSoft; }
+		}
 		sigmaT = G*epsT; /* trial stress */
-		Real yieldSigmaT = max((Real)0.,undamagedCohesion*(1-omega)-sigmaN*tanFrictionAngle); /* Mohr-Coulomb law with damage */
+		Real yieldSigmaT = std::max((Real)0.,undamagedCohesion*(1-omega)-sigmaN*tanFrictionAngle); /* Mohr-Coulomb law with damage */
 		if (sigmaT.squaredNorm() > yieldSigmaT*yieldSigmaT) {
 			Real scale = yieldSigmaT/sigmaT.norm();
 			sigmaT *= scale; /* stress return */
@@ -588,7 +596,7 @@
 		if (!phys->isCohesive) continue;
 		bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum += (1-phys->relResidualStrength); // bodyStats[id1].epsPlSum += phys->epsPlSum;
 		bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum += (1-phys->relResidualStrength); // bodyStats[id2].epsPlSum += phys->epsPlSum;
-		maxOmega = max(maxOmega,phys->omega);
+		maxOmega = std::max(maxOmega,phys->omega);
 		avgRelResidual += phys->relResidualStrength;
 		nAvgRelResidual += 1;
 		for (int i=0; i<3; i++) {

=== modified file 'pkg/dem/ConcretePM.hpp'
--- pkg/dem/ConcretePM.hpp	2016-04-23 17:54:37 +0000
+++ pkg/dem/ConcretePM.hpp	2016-05-13 13:07:22 +0000
@@ -104,6 +104,7 @@
 		((bool,neverDamage,false,,"If true, no damage will occur (for testing only)."))
 		((Real,epsCrackOnset,NaN,,"Limit elastic strain [-]"))
 		((Real,relDuctility,NaN,,"relative ductility of bonds in normal direction"))
+		((Real,equivStrainShearContrib,0,,"Coefficient of shear contribution to equivalent strain"))
 		((int,damLaw,1,,"Law for damage evolution in uniaxial tension. 0 for linear stress-strain softening branch, 1 (default) for exponential damage evolution law"))
 		((Real,dmgTau,((void)"deactivated if negative",-1),,"Characteristic time for normal viscosity. [s]"))
 		((Real,dmgRateExp,0,,"Exponent for normal viscosity function. [-]"))
@@ -168,6 +169,7 @@
 			((Real,epsCrackOnset,NaN,,"strain at which the material starts to behave non-linearly"))
 			((Real,relDuctility,NaN,,"Relative ductility of bonds in normal direction"))
 			((Real,epsFracture,NaN,,"strain at which the bond is fully broken [-]"))
+			((Real,equivStrainShearContrib,NaN,,"Coefficient of shear contribution to equivalent strain"))
 			((Real,dmgTau,-1,,"characteristic time for damage (if non-positive, the law without rate-dependence is used)"))
 			((Real,dmgRateExp,0,,"exponent in the rate-dependent damage evolution"))
 			((Real,dmgStrain,0,,"damage strain (at previous or current step)"))