cpm model tensile contact breakage

Asked by chanaka Udaya on 2019-10-30

Hi all,

The problem is regarding CPM model when it subjected to tension failure.

This contact law has exponential softening damage in tension. Therefore,after some time, particles are subjected to fully damage during the simulation.at this point, value of the damage variable is equal to 1.

My problem is ,what will happen to the tensile contact after when it is subjected to fully damage(omega=1)?

Referring to source code ConcretePM.hpp, I can see following line,(line 281)

((Real,omegaThreshold,((void)">=1. to deactivate, i.e. never delete any contacts",1.),,"damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞"))

and ConcretePM.cpp (line 420-429)

/* handle broken contacts */
 if (epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)) {
  /* if (isCohesive) { */
    const shared_ptr<Body>& body1 = Body::byId(I->getId1(),scene), body2 = Body::byId(I->getId2(),scene); assert(body1); assert(body2);
    const shared_ptr<CpmState>& st1 = YADE_PTR_CAST<CpmState>(body1->state), st2 = YADE_PTR_CAST<CpmState>(body2->state);
   /* nice article about openMP::critical vs. scoped locks: http://www.thinkingparallel.com/2006/08/21/scoped-locking-vs-critical-in-openmp-a-personal-shootout/ */
   { boost::mutex::scoped_lock lock(st1->updateMutex); st1->numBrokenCohesive += 1; /* st1->epsPlBroken += epsPlSum; */ }
   { boost::mutex::scoped_lock lock(st2->updateMutex); st2->numBrokenCohesive += 1; /* st2->epsPlBroken += epsPlSum; */ }
  /* } */
  return false;
 }

I couldn't understand the lines in .cpp file with assert(body1), assert(body2).since if we want to use assert function in c++, there should be an expression.

Thank you
Chanaka

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-10-31
Last query:
2019-10-31
Last reply:
2019-10-31
Jan Stránský (honzik) said : #1

Hello,

> ... exponential softening ... value of the damage variable is equal to 1 ...

because of the exponential law, value 1 is the limit, but the value is never 1 (although could be very close to 1).
Maybe for very high tensile strain because of round-off errors it can be 1 or slightly higher, but it should not be usual case.

> My problem is ,what will happen to the tensile contact after when it is subjected to fully damage(omega=1)?

nothing (see above).
If the damage is higher than omegaThreshold, the contact is deleted.
(the default value omegaThreshold=1 means that contacts are not deleted even for very high damage - because it should not be 1 or higher)

> I couldn't understand the lines in .cpp file with assert(body1), assert(body2).since if we want to use assert function in c++, there should be an expression.

body1 is a pointer to a Body instance. It can be nullptr, the assertion is that body1 and body2 are not nullptr.

cheers
Jan

chanaka Udaya (chanaka-udaya) said : #2

Hello Jan,

many thanks for the detailed answer.

If I want to delete the contact after damage is reaching to 0.9999, Is it okay to change the source code like below ?,

 ConcretePM.hpp,(line 281)

((Real,omegaThreshold,((void)">=1. to deactivate, i.e. never delete any contacts",0.9999),,"damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞"))

I have changed the omegaThreshould to 0.9999 and in .cpp file I did not change anything.

 if (epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)) {
  /* if (isCohesive) { */
    const shared_ptr<Body>& body1 = Body::byId(I->getId1(),scene), body2 = Body::byId(I->getId2(),scene); assert(body1); assert(body2);
    const shared_ptr<CpmState>& st1 = YADE_PTR_CAST<CpmState>(body1->state), st2 = YADE_PTR_CAST<CpmState>(body2->state);
   /* nice article about openMP::critical vs. scoped locks: http://www.thinkingparallel.com/2006/08/21/scoped-locking-vs-critical-in-openmp-a-personal-shootout/ */
   { boost::mutex::scoped_lock lock(st1->updateMutex); st1->numBrokenCohesive += 1; /* st1->epsPlBroken += epsPlSum; */ }
   { boost::mutex::scoped_lock lock(st2->updateMutex); st2->numBrokenCohesive += 1; /* st2->epsPlBroken += epsPlSum; */ }
  /* } */
  return false;
 }

After compiling, I got a following message for tension between two particles after some time,

AttributeError: 'NoneType' object has no attribute 'normalForce'

That means the contact is deleted?

Jan Stránský (honzik) said : #3

Hi,

> If I want to delete the contact after damage is reaching to 0.9999, Is it okay to change the source code like below ?,

that is one option.
Far easier option is not to touch the source code at all and specify the value in your python script :-)
###
O.engines = [
   ...
      Law2_ScGeom_CpmPhys_Cpm(omegaThreshold=whatever),
   ...
]
###

> That means the contact is deleted?

maybe. Please provide thescript you are using.

cheers
Jan

chanaka Udaya (chanaka-udaya) said : #4

Hi Jan,

I have changed the line below in source code and after compiling ran the following script.

ConcretePM.hpp,(line 281)

((Real,omegaThreshold,((void)">=1. to deactivate, i.e. never delete any contacts",0.9999),,"damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞"))

here is the my script.

from yade import plot,qt
nh=O.materials.append(CpmMat(young=20e9,frictionAngle=atan(0.5),poisson=.5,density=2600,sigmaT=1.5e6,relDuctility=1.5,epsCrackOnset=5e-5,isoPrestress=0))
s1=utils.sphere((0,0,0),radius=1,color=[1,0,0],material=nh)
s2=utils.sphere((0,0,2),radius=1,color=[0,1,0],material=nh)

O.bodies.append(s1)
O.bodies.append(s2)

Gl1_Sphere.quality=3
# engines definition
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.3,label='bo1s')]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.3,label='ig2ss')],
        [Ip2_CpmMat_CpmMat_CpmPhys()],
        [Law2_ScGeom_CpmPhys_Cpm()]
    ),
    TranslationEngine(ids=[s2.id],translationAxis=[0,0,1],velocity=0.01*100), #shear+tension
    TranslationEngine(ids=[s1.id],translationAxis=[0,0,-1],velocity=0.01*100),
    NewtonIntegrator(damping=0.4),
    PyRunner(iterPeriod=1,command='addPlotData()',label='plotData'),
]
O.dt=0.004*PWaveTimeStep()
O.step()
bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1.
def addPlotData():
    plot.addData(i=O.iter,
                 displacement_n=O.bodies[1].state.refPos[2]-O.bodies[1].state.pos[2],
                 NF=O.interactions[0,1].phys.normalForce.norm(),
                 SF=O.interactions[0,1].phys.shearForce.norm(),
                 )

plot.plots={'i':('NF'),
           'displacement_n ':('NF'),
    'i ':('SF'),
   }
plot.plot()
O.run()

Thanks

Best Jan Stránský (honzik) said : #5

According to the code and error, the interaction is deleted (O.interactions[0,1] returned None)
cheers
Jan

chanaka Udaya (chanaka-udaya) said : #6

Thanks Jan Stránský, that solved my question.