Is possible to plot the unbreakable bonds during UTS in paraview

Asked by kalogeropoulos on 2021-01-20

Hello to everybody, i hope everyone is fine during Convid!!

To be honest, i beleive the title of my question says everything.
I am conducting UTS simulations and i want to process the reuslts in paraview.
I dont know how to plot the unbreakable zone of bonds in paraview?
Is any way fir something like that to be achieved?

Thank you.
I am showing the code below...

# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-
from yade import ymport, plot

#### Simulation of a uniaxial test (compression or tension depending on the sign of the loading rate "rate" defined below)
# everything is documented at https://yade-dem.org/doc/index-toctree.html -> use the quick serach bar on the left of the page

################# SIMULATIONS DEFINED HERE

#### Simulation Control
rate=0.01 # deformation rate (0.003 for tension, -0.03 for compression)
iterMax=10000 # maximum number of iterations
saveVTK=4000 # saving output files for paraview
OUT='tensionTest_r0.01' # output files

intR=1.00 # allows near neighbour interaction (can be adjusted for every packing)
DENS=5265 # 2640 toy petrwmatos# could be adapted to match material density: dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> packing porosity as to be computed?
YOUNG=6e2
FRICT=13
ALPHA=2.5
TENS=2e6
COH=1.5*TENS

#### material definition
def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson=ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
D=5e-2
H = 2.5*D

pred=pack.inCylinder((0,0,0),(0,H,0),D)

sps=SpherePack()
sp=pack.randomDensePack(pred,spheresInCell=10000,radius=1e-3,rRelFuzz=.34,memoizeDb='/tmp/test.sqlite',returnSpherePack=True)
sp.toSimulation(material=sphereMat) # up to you to define another sample here, e.g., with randomDensePack or anything else.

#### Particle Characteristics
R=0
Rmax=0
nbSpheres=0.
for o in O.bodies:
 if isinstance(o.shape,Sphere):
   nbSpheres+=1
   R+=o.shape.radius
   if o.shape.radius>Rmax:
     Rmax=o.shape.radius
Rmean=R/nbSpheres
print "\n"
print ('nbSpheres=',nbSpheres,' | Rmean=',Rmean,'| Rmax=',Rmax,)
print "\n"

#### Set up boundary conditions (see https://yade-dem.org/doc/yade.utils.html?highlight=uniaxialtestfeatures#yade.utils.uniaxialTestFeatures)
bb=utils.uniaxialTestFeatures()
negIds,posIds,longerAxis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']

################# ENGINES DEFINED HERE

O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom')],
  [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
  [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
 ),
 UniaxialStrainer(strainRate=rate,axis=longerAxis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=1,blockRotations=1,setSpeeds=0,stopStrain=0.1,dead=1,label='strainer'),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
 NewtonIntegrator(damping=0.4,label='newton'),
 PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
 PyRunner(iterPeriod=int(100), command='Stop()',label='data'),
        VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
]

################# RECORDER DEFINED HERE

def recorder():
    yade.plot.addData({'i':O.iter,
         'eps':strainer.strain,
         'sigma':strainer.avgStress,
         'tc':interactionLaw.nbTensCracks,
         'sc':interactionLaw.nbShearCracks,
         'te':interactionLaw.totalTensCracksE,
         'se':interactionLaw.totalShearCracksE,
         'unbF':utils.unbalancedForce()})
    plot.saveDataTxt(OUT)

def Stop():
    sigma=plot.data['sigma']
    if abs(strainer.avgStress) <= 0.75*max(sigma):
        O.pause()

# if you want to plot during simulation
plot.plots={'eps':('sigma')}
plot.plot()

################# PREPROCESSING
print " Number of particles: %.2f" %(len(O.bodies))
print " Porosity of sample is: %.4f" %(utils.porosity())
print "\n"

#### manage interaction detection factor during the first timestep and then set default interaction range ((cf. A DEM model for soft and hard rock, Scholtes & Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

# reinforce bonds at boundaries for tensile tests
dim=utils.aabbExtrema()
if rate>0:
 layerSize=0.15
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   if ( o.state.pos[longerAxis]<(dim[0][longerAxis]+layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) or ( o.state.pos[longerAxis]>(dim[1][longerAxis]-layerSize*(dim[1][longerAxis]-dim[0][longerAxis])) ) :
    o.shape.color=(1,1,1)
 for i in O.interactions:
  if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
   if O.bodies[i.id1].shape.color==(1,1,1) or O.bodies[i.id2].shape.color==(1,1,1) :
    i.phys.FnMax*=100
    i.phys.FsMax*=100

#### coordination number after the creation of contacts at first step. (Before = 0)
print " Number of Interactions per Particle: %.2f" %(utils.avgNumInteractions())#CoordinationNumber
print "\n"
#### coordination number verification
numSSlinks=0
numCohesivelinks=0
for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      numSSlinks+=1
    if i.phys.isCohesive :
      numCohesivelinks+=1

print ("nblinks=", numSSlinks, " | nbCohesivelinks=", numCohesivelinks, " || Kcohesive=", 2.0*numCohesivelinks/nbSpheres)
print "\n"
################# SIMULATION REALLY STARTS HERE
strainer.dead=0
#O.run(iterMax)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Luc Scholtès
Solved:
2021-01-21
Last query:
2021-01-21
Last reply:
2021-01-21
Luc Scholtès (luc) said : #1

Hi,

The 'jcfpm' key in the VTKRecorder should produce '.intr.vtp' files that you can then open with Paraview. If that's not the case (it should not), you need to add the 'intr' key in the recorders' list.

Once you have these '.intr.vtp' files, you can open them in Paraview and then select the 'isCohesive' option to see the "unbroken" bonds (the cohesive ones).

Let us know if that's what you were looking for.

Luc

kalogeropoulos (antoniskal) said : #2

Well I didn't exactly mean that.

I was telling about the fixed layer of spheres during UTS.
Is it possible somehow to have for example, with black colour the fixed spheres and in yellow the rest sample during the UTS in paraview?

Best Luc Scholtès (luc) said : #3

Ah.. sorry I misread your question.

For the bonds, I don't see anything that could help but it is absolutely feasible for the spheres. Just color the spheres differently in the simulation (e.g., spheres in layer in red, color=(1,0,0) and spheres outside in blue, color= (0,1,0)) and then add the key 'colors" in the VTKRecorder. You should be able to play with the color of the spheres when opening the .sphere.vtu in Paraview. Actually, given the current settings (color of spheres in your simulation), you should already be able to distinguish the layer in Paraview (just add 'colors' in the VTKRecorder).

Luc

kalogeropoulos (antoniskal) said : #4

Thanks Luc Scholtès, that solved my question.