Force in a facet

Asked by kalogeropoulos on 2018-10-17

Hello launchpad,

I am new to Yade. I have read the doc, and the examples.
I am trying to simulate a rock cutting and record :
1) the forces acting on the sample, and
2) the cracks

Here is my code:

#!/usr/bin/python
# -*- coding: utf-8 -*-
from __future__ import division

from yade import plot,pack,timing, utils, geom
import time, sys, os, copy

#import matplotlib
#matplotlib.rc('text',usetex=True)
#matplotlib.rc('text.latex',preamble=r'\usepackage{concrete}\usepackage{euler}')

# default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
 young=24e9, #kn
 poisson=.2, #ks

 sigmaT=3.5e6,
 frictionAngle=atan(0.8),
 epsCrackOnset=1e-4,
 relDuctility=30,

 intRadius=1.5,
 dtSafety=.8,
 damping=0.4,
 #strainRateTension=.05,
 #strainRateCompression=.5,
 #setSpeeds=True,
 # 1=tension, 2=compression (ANDed; 3=both) den xreiazontai, prokeitai gia dokimh kophs
 #doModes=2,

 specimenLength=1,
 sphereRadius=0.99e-3, #0.35e-3

 # isotropic confinement (should be negative)
 isoPrestress=0,
)

from yade.params.table import *

if 'description' in O.tags.keys(): O.tags['id']=O.tags['id']+O.tags['description']
#material properties
concreteId=O.materials.append(CpmMat(young=young,frictionAngle=frictionAngle,poisson=poisson,density=4800,sigmaT=sigmaT,
     relDuctility=relDuctility,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress))

sps=SpherePack()
sp=pack.randomDensePack(pack.inAlignedBox((-75e-3*specimenLength,-27e-3*specimenLength,-6.35e-3*specimenLength),(75e-3*specimenLength,27e-3*specimenLength,6.35e-3*specimenLength)),spheresInCell=2000,radius=sphereRadius,memoizeDb='/tmp/triaxPackCache.sqlite',returnSpherePack=True)
sp.toSimulation(material=concreteId)

bx = geom.facetBox(center=(-80e-3,27e-3,0.0),extents=(3e-3,-5e-3,6.35e-3),orientation=Quaternion((0,0,1),-pi/36), wallMask=(2+8))
O.bodies.append(bx)
for facet in bx:
 facet.state.blockedDOFs='xyzXYZ'
 facet.state.vel=(20.0,0,0)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Facet_Aabb()],verletDist=.05*sphereRadius),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'), Ig2_Facet_Sphere_ScGeom()],
  [Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)],
  [Law2_ScGeom_CpmPhys_Cpm()],
 ),
 NewtonIntegrator(damping=damping,label='damper'),
 PyRunner(command='addPlotData()',realPeriod=2), #CpmStateUpdater(realPeriod=.5),
]

def addPlotData():
   for facet in bx:
      #Fn = facet.shape.normal
      Fn=O.forces.f(facet.id)[2]
      disp = facet.state.pos[2]
   plot.addData(i=O.iter,Fn=Fn,disp=disp)
plot.plots={'displacement':('Fn',)}
plot.plot()

O.dt=dtSafety*PWaveTimeStep()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
kalogeropoulos
Solved:
2018-10-18
Last query:
2018-10-18
Last reply:
2018-10-17

This question was reopened

Hi, Welcome to Yade.
Do you have a question?
Bruno

kalogeropoulos (antoniskal) said : #2

I am sorry, i forgot to put the question.

How can i plot the forces that acting on the facets?

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

Why was it solved?

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

Hello

> I am new to Yade

welcome :-)

> How can i plot the forces that acting on the facets?

In your current approach in addPlotData, you define Fn and dspl for all facets, but after the cycle, only the last Fn is saved to plot.data. To use total force, either sum it in the for cycle
#
Fn = 0.0
for facet in facets:
   Fn += O.forces.f(facet.id)[2]
#
or use python shorthand: Fn = sum(O.forces.f(facet.id)[2] for facet in facets)

you then plot is either in Yade or you can plot.saveDataTxt and plot it externally (using excel, gnuplot,...)

Also note that you prescribe velocity of facets in x direction but measure z force, is it ok?

cheers
Jan

kalogeropoulos (antoniskal) said : #6

Thnaks a lot.

I did it as Jan Stránský (honzik) mentioned and worked. Also i change the Fn = sum(O.forces.f(facet.id)[2] to
Fn = sum(O.forces.f(facet.id)[0] to record the Force in x-direction.

Thanks a lot.