Yade

I got a wrong damage ratio

Asked by Huihuang Xia

Hi,

https://answers.launchpad.net/yade/+question/432617, in this question, damage ratio was defined and implemented in YADE.
According to this definition, I built a impact crushing simulation using damage ratio in YADE, but my code got a wrong result, as damage ratio in this simulation was greater than 1.

Here is my code:

#!/usr/bin/python
# -*- coding: utf-8 -*-

#Created by Huihuang Xia from Huaqiao University,Xiamen,P.R.C
#My E-Mail:<email address hidden>

###########################
# IMPORT MODULES
###########################

from yade import pack
from yade import plot
from yade import qt

###########################
# DEFINE MATERIALS
###########################

rock=CohFrictMat(young=5.98e7,poisson=0.3,alphaKr=3000,alphaKtw=3000,density=2678,frictionAngle=0.5,isCohesive=True,normalCohesion=7.9e6,shearCohesion=7.9e6,momentRotationLaw=True)
O.materials.append(rock)
steel=CohFrictMat(young=3.06e11,poisson=0.29,density=7861,frictionAngle=0.545,normalCohesion=0,shearCohesion=0)
O.materials.append(steel)

#################################
# CREATE SAMPLE & RIGID_WALL
#################################
wall=O.bodies.append(geom.facetBox(center=(0,0,0),extents=(0.015,0.015,0.0005),color=(1,1,0),material=steel,fixed=True))
pred=pack.inSphere(center=(0,0,0.006),radius=0.005)
assembly=pack.randomDensePack(pred,radius=0.0002,rRelFuzz=0.5,spheresInCell=2500,material=rock)
O.bodies.append(assembly)
for b in assembly:
 b.state.vel=(0,0,-20)

###########################
# ENGINES
###########################
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.0),Bo1_Facet_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=1.0),Ig2_Facet_Sphere_ScGeom6D()],
        [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True)],
        [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
    VTKRecorder(fileName='post/impact-',recorders=['all'],iterPeriod=250),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=5,timestepSafetyCoefficient=0.8,defaultDt=PWaveTimeStep()),
    NewtonIntegrator(gravity=(0,0,-9.81)),
    PyRunner(command='damageRatio()',realPeriod=5),
    PyRunner(command='addPlotData()',realPeriod=5),
]
O.trackEnergy=True
O.step()

###########################
# DEFINE FUNCATIONS
###########################

#define a function to save damage ratio
global sumCohBonds
sumCohBonds=0
for b in O.interactions:
    if b.phys.cohesionBroken==True:
        continue
    sumCohBonds+=1 # calculate total intact bonds

numBroCohBonds=0
def damageRatio():
    global numBroCohBonds
    for br in O.interactions:
        if br.phys.cohesionBroken==False:
         continue
        numBroCohBonds+=1 # calculate broken bonds
    damageRatio=numBroCohBonds/float(sumCohBonds)
    print sumCohBonds,numBroCohBonds
    plot.saveDataTxt('data/damageratio.txt')
def addPlotData():
    damageRatio=numBroCohBonds/float(sumCohBonds)
    plot.addData(t=O.time,dr=damageRatio)
plot.plots={'t':('dr')}
plot.plot()

#set an optimal timestep
O.dt=utils.PWaveTimeStep()
O.usesTimeStepper=True

#3D view and controller
qt.View()
qt.Controller()

and here is the damage ratio vs time:

# dr t
0.0 9.49679622677e-09
0.0044209618287 1.60495856089e-06
0.0155272805693 4.9193404254e-06
0.0266335993099 8.39516779415e-06
0.0362303213284 1.10352770982e-05
0.0474444684063 1.44826140458e-05
0.0585507871469 1.73886336283e-05
0.0679318524908 2.02186788283e-05
0.0788225145568 2.3419099093e-05
0.091654086694 2.71228495603e-05
0.104485658831 3.03992441997e-05
0.117856372655 3.34097285591e-05
0.130687944792 3.48247511657e-05
0.144058658615 3.669561998e-05
0.158291999137 3.75218412351e-05
0.172633167997 3.86994439409e-05
0.190748328661 4.02569184773e-05
0.210696571059 4.18618769955e-05
0.236144058659 4.52902203506e-05
0.278089281863 4.98581792485e-05
0.339120120768 5.47350982335e-05
0.42009920207 5.83200136999e-05
0.520379555747 6.13408059573e-05
0.640931636834 6.41705226598e-05
0.776148371792 6.70709729927e-05
0.946409316368 7.14288492846e-05
1.15074401553 7.49846615977e-05
1.39076989433 7.87721800155e-05
1.67392710804 8.25338810098e-05
1.99611817986 8.60126620602e-05
2.32769031702 8.71703956744e-05
2.67360362303 8.81834109127e-05
3.04539573 9.07206891265e-05
3.46236791029 9.43729049081e-05
3.92160879879 9.80223033349e-05
4.41718783696 0.000101386989899
4.94899719646 0.000104757207048
5.50323485012 0.000107345725132
6.08561569981 0.000110695198656
6.67565236144 0.000111623153105
7.272697865 0.000112554955106
7.87524261376 0.000113345795924
8.48781539789 0.000116016886332
9.10739702394 0.000118131382236
9.72881173172 0.000120064040351
10.3572352814 0.000121999979727
10.989756308 0.000123849996843
11.6188268277 0.000126003034127

PS: this simulation got no error and my Yade (V.2016.06a)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Luc Scholtès
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said : 		#1

Hello,

It looks like you are summing through broken bonds at each timestep and adding it to the sum from the previous step. If you want your damage ratio to be relative to each timestep, the number of broken bonds needs to be reset before each timestep.

So I would put the numBroCohBonds=0 inside the damageRatio() function. This way,the summation is reset each time damageRatio() is called.

Robert

Revision history for this message
Huihuang Xia (huihuangxia) said : 		#2

Hi Robert,

Thanks for your answer, but put the numBroCohBonds=0 inside the damageRatio() function is reallly a big mistake. You can refer to https://answers.launchpad.net/yade/+question/432617, in this question, it was stated that put the numBroCohBonds=0 inside the damageRatio() function led to a initialization of numBroCohBonds to zero each time.
Damage ratio is the fraction of bonds broken relative to the total number of bonds originally created. In each time step, some new bonds broke, and number of bonds broken in damage ratio must be a sum of broken bonds in all time step before current time step. Thus, I still need an answer.

Huihuang Xia

Revision history for this message
Robert Caulk (rcaulk) said : 		#3

Hello Huihuang,

I recommend you read the final post of that thread more closely ;-). You [i]want[/i] to initialize numBroCohBonds to zero each time, because each time this function is called you are cycling through all of the interactions and summing all the brokenbonds.

Maybe you misunderstood what I was recommending. Here it is:

def damageRatio():
    global numBroCohBonds
    numBroCohBonds=0
    for br in O.interactions:
        if br.phys.cohesionBroken==False:
         continue
        numBroCohBonds+=1 # calculate broken bonds
    damageRatio=numBroCohBonds/float(sumCohBonds)
    print sumCohBonds,numBroCohBonds
    plot.saveDataTxt('data/damageratio.txt')

Best,

Robert

Revision history for this message
Huihuang Xia (huihuangxia) said : 		#4

Hi Robert,

Thanks for your patience, but your method also got a wrong result, as damage ratio was greater than 1 and number of broken bonds could not increase gradually, which violated the definition of damage ratio.

Revision history for this message
Luc Scholtès (luc) said : 		#5

I see 2 potential problematic reasons why you get your damage ration >1:

1) you don't use the same condition for checking the nb of cohesive contacts at initialization and during the simulation -> Why?

initialization:
if b.phys.cohesionBroken==True:
continue
sumCohBonds+=1

during simulation:
if br.phys.cohesionBroken==False:
         continue
numBroCohBonds+=1

Revision history for this message
Best Luc Scholtès (luc) said : 		#6

Sorry, for wrong use of keyboard shortcut... I try another time here:

I see 2 potential reasons why you get your damage ratio >1:

1) you don't use the same condition for checking the nb of cohesive contacts at initialization and during the simulation:

initialization:
if cohesionBroken==True:
   continue
sumCohBonds+=1

during simulation:
if cohesionBroken==False:
   continue
numBroCohBonds+=1

Why don't you use the same condition for both cases?

2) the test you use during the simulation might cause problem because you sum over interactions and look for those with cohesionBroken==True and this paramater is set to true by default, even for frictional interactions (if I understood correctly). Doing so, the newly created interactions are counted as being broken despite the fact that they appeared as a result of the interaction of particles after breakage of you specimen.

For me, the solution would be to use the same condition for each case:

if cohesionBroken==False:
   numCohBonds+=1

and to define InitialNumCohBonds in the initial loop and currentNumCohBonds in the current loop so that the damage ratio would be:

damageRatio=currentNumCohBonds/initialNumCohBonds

If, after that, you still have a damage ratio > 1, there might be a problem somewhere else... For instance, I have never used this contact law before but are you sure about the way you use the setCohesionNow flag? In your case, you define it equal to true in the engine list -> Does it mean that it will be always true or just for the first iteration of your simulation? From my understanding, the way you define it makes it true for all the duration of the simulation hence, you create cohesive contacts every time particles come into contact hence, your damage ratio > 1.

Revision history for this message
De zhang (dzlearnyade) said : 		#7

change the code as flowing:
def addPlotData():
     damageRatio=numBroCohBonds/float(numBroCohBonds+sumCohBonds)
     plot.addData(t=O.time,dr=damageRatio)
haha~~

Revision history for this message
Huihuang Xia (huihuangxia) said : 		#8

Thanks Luc Scholtès, that solved my question.