Hello,
I am a new YADE user, trying to calibrate my DEM model
by performing uniaxial compression tests with UniaxialStrainer.
I need to compute transverse strain in a cylindrical specimen,
in order to match it with macroscopic poisson’s ratio.
Any ideas on how i could do it?
I think its not “part” of uniaxial strainer.
Thanks

Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2019-11-28
Last query:
2019-11-28
2019-11-28
 Jérôme Duriez (jduriez) said on 2019-11-14: #1

Hi,

No, it's not part of UniaxialStrainer which considers only axial quantities. You need to measure yourself the transversal dimensions of your specimen. First possibility could be to use aabbDim()

 Jan Stránský (honzik) said on 2019-11-14: #2

Hello,

> I think its not “part” of uniaxial strainer.

exactly. Because computation of it is not very straightforward. E.g. you can average "particle transverse strain" as newDistanceFromAxis/originalDistanceFromAxis-1. Then particles near the axis would be problematic because their original length is low (or possibly 0) and would corrupt the results.
Also it may happen that boundary particles "fly away", also corrupting the result.
So you can e.g. somehow choose a "central ring" and do computation on these particles.

Apart from computing transverse strain from uniaxial stress test, there is another option how to compute Poisson's ratio.
Poisson's ratio can also be computed by comparing Young's modulus (uniaxial stress modulus) to oedometric (uniaxial strain) modulus and/or bulk ("triaxial") modulus. See related discussion [1].

You can use both methods (transverse strain and another loading) and compare the results of Poisson's ratio :-)

cheers
Jan

 AnnaMaria (annm5) said on 2019-11-14: #3

Thanks for your response Jérôme Duriez (jduriez).
I read in [1], as Jan Stránský (honzik) suggested that you have
to to define which particles are "suitable" for the computation.
That very near the axis had negative influence, as well as those "crushed away" from the surface and flying away. I choose those roughly half-way from the axis to surface.
Can I do this with aabbDim()? Choose those half-away?
Access a region within the cylinder?
utils.aabbDim(cutoff=0.0, centers=False)

Thanks

 AnnaMaria (annm5) said on 2019-11-14: #4

Thanks Jan Stránský (honzik), for your response
it was while i was typing a comment on what Jérôme Duriez (jduriez) said
and i did not see it in time.

Im not sure how i can access this ring. With aabbDim?

Thanks

 Jan Stránský (honzik) said on 2019-11-14: #5

> Can I do this with aabbDim()?

no, aabbDim computes aabb (axis aligned bounding box) dimensions of the sample (roughly the "maximum dimensions along axes").

To choose particles within a "ring" defined by radii r1 and r2, do:
### (not tested)
halfwayParticles = []
for b in O.bodies:
x,y,z = b.state.pos # maybe b.state.refPos
d = sqrt(pow(x,2)+pow(y,2)) # assuming z axis
if d > r1 and d < r2:
halfwayParticles.append(b)
###

cheers
Jan

 AnnaMaria (annm5) said on 2019-11-14: #6

Thanks Jan Stránský (honzik)
i will try what you said ;)

 Bruno Chareyre (bruno-chareyre) said on 2019-11-18: #7

An alternative approach, which may fit your needs or not, would be to use TriaxialStressController by prescribing a very small confinement on the other two directions. That would prevent particles from flying away and it would keep the boundaries in contact with the specimen, which makes defining transverse strain straightforward.
If the confinement is small enough it should approach the uniaxial response.
Cheers
Bruno

 Jan Stránský (honzik) said on 2019-11-18: #8

Bruno's answer would work for prismatic shape. If you have cylinder, it would not work.

@Anna: what is the shape fo your specimen?

cheers
Jan

 AnnaMaria (annm5) said on 2019-11-19: #9

My specimen is cylindrical.
I am still trying to make it work. (the way that Jan Stránský (honzik) suggested)

 AnnaMaria (annm5) said on 2019-11-28: #10

Hello everyone again:
I have created the hollow cylinder as Jan Stransky suggested, (MWE):
for p in O.bodies:
x=p.state.refPos[0]
z=p.state.refPos[2]
d=sqrt(pow(x,2)+pow(z,2))
if d>r1 and d<r2:
p.shape.color=(1,1,1)
#Mesure state.refPos of the particles, within this domain.
xnew=p.state.refPos[0]
znew=p.state.refPos[2]
dnew=sqrt(pow(xnew,2)+pow(znew,2))
Then I am trying to record the lateral strain:
def recorder():
global d
global dnew
elat=((d)-(dnew))/(d)
but it does not work. I am not getting any values!
I can’t understand what i have been doing wrong all this time.
Could you help me. Thanks in advance.

 Jan Stránský (honzik) said on 2019-11-28: #11

Please provide a real MWE (W=working, a complete script which we can test ourselves).
Thanks
Jan

 AnnaMaria (annm5) said on 2019-11-28: #12

Thanks Jan,
This is a working example.
I want to plot in the positive domain of the Cartesian coordinate system the stress vs strain, and in the negative the lateral strain. But right now it does not seem to work, even if only plot lateral strain. Stress seems to plot ok.

#### Simulation Control
rate=-1 #deformation rate

#### Material
# default parameters or from table
Specimen_Height= 0.13,
intR=1.5,
DENSITY=3000,
YOUNG=60e9,
FRICTION_ANGLE=50,
POISSON=1,
TENS=5e6,
COHESIO=40e6,
r1=0.010,
r2=0.020,
h1=0.061,
h2=0.061,
)

#### material definition
Sample = O.materials.append(JCFpmMat(
young=YOUNG,
poisson=POISSON,
cohesion=COHESIO,
tensileStrength=TENS,
density = DENSITY,
label='spheres'
))
#### Specimen Creation
sp=pack.randomDensePack(
spheresInCell=1000,
material=Sample,
rRelFuzz=0.12,
memoizeDb='/tmp/Cpack.db',
returnSpherePack=True,
)
sp.toSimulation()
for p in O.bodies:
x=p.state.refPos[0]
z=p.state.refPos[2]
d=sqrt(pow(x,2)+pow(z,2))
if d>r1 and d<r2:
p.shape.color=(1,1,1)
xnew=p.state.refPos[0]
znew=p.state.refPos[2]
dnew=sqrt(pow(xnew,2)+pow(znew,2))

#### boundary condition (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,label='interactionLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.5, defaultDt=utils.PWaveTimeStep()),
NewtonIntegrator(damping=0.2,label='newton'),
PyRunner(iterPeriod=int(100),initRun=True,command='recorder()',label='data'),
]

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

def recorder():
global d
global dnew
elat=(d-dnew)/d
'eps':abs(strainer.strain),
'elat':elat,
'sigma':abs(strainer.avgStress),})

#### Plot During Simulation
plot.plots={'elat':('sigma')}
plot.plot()

O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

 Jan Stránský (honzik) said on 2019-11-28: #13

1) I have no problem running the script, plot.data shows the saved data.

2) Python/Yade reads the script and executes it line by line. In your case, it loops over the bodies and computes the d and dnew ONLY once at the beginning, then the value is unchanged during the whole simulation.
Put the code inside a function and call that function inside recorder() before plot.addData()

3) furthermore, d and dnew are just values form the last single body.
You want probably some average (maybe wighted?) of the values.

4) expressions for d and dnew does not differ, probably dnew should be computed using state.pos instead of state.refPos?

cheers
Jan

 AnnaMaria (annm5) said on 2019-11-28: #14

Hi Jan,
I think I am getting close
2) I Changed the script as follows:
After creating the packing:
#sp.toSimulation()
def lateral():
global elat
global d
global dnew
for p in O.bodies:
x=p.state.refPos[0]
z=p.state.refPos[2]
d=sqrt(pow(x,2)+pow(z,2))
if d>r1 and d<r2:
p.shape.color=(1,1,1)
xnew=p.state.pos[0]
znew=p.state.pos[2]
dnew=sqrt(pow(xnew,2)+pow(znew,2))
elat=(d-dnew)/d
return elat

and then:
def recorder():
lateral()

as you said.
The thing is that the function above does not seem to work correct. The spheres have the same color, while before they were white (color=(1,1,1)). Maybe the distance measurement is not working as well.
Also what do you mean:
3) furthermore, d and dnew are just values form the last single body.
You want probably some average (maybe wighted?) of the values.
Like this: elat=(avg(d)-avg(dnew))/avg(d)

 Jan Stránský (honzik) said on 2019-11-28: #15

5)
Use more than one space for indentation. Here your recorder would be
###
def recorder():
lateral()
###
meaning plot.addData outside recorder(), which is probably not what you want

6)
> The spheres have the same color
no, there is one sphere with different color.
In laterla(), you return directly when the very first suitable body is found.

6)+3) solution
###
def lateral():
elatTot = 0.0 # sum of elat of all suitable particles
nTot = 0 # number of suitable particles
for p in O.bodies:
x=...; z=...; d=...
if d>r1 and d<r2:
xnew=...; znew=...; dnew=...
elat=(d-dnew)/d
elatTot += elat
nTot += 1
return elatTot / nTot # average elat
###

cheers
Jan

 AnnaMaria (annm5) said on 2019-11-28: #16

Jan Thank you very much for your effort!!!!
I really appreciate it.
I have only changed elatTot / nTot to elat_avg=elatTot / nTot
And then return elat_avg # average elat so that I can plot.

def lateral():
elatTot = 0.0 # sum of elat of all suitable particles
nTot = 0 # number of suitable particles
for p in O.bodies:
x=p.state.refPos[0]
z=p.state.refPos[2]
d=sqrt(pow(x,2)+pow(z,2))
if d>r1 and d<r2:
xnew=p.state.pos[0]
znew=p.state.pos[2]
dnew=sqrt(pow(xnew,2)+pow(znew,2))
elat=(d-dnew)/d
elatTot += elat
nTot += 1
elat_avg=elatTot / nTot
return elat_avg # average elat

def recorder():
lateral()
plot.saveDataTxt(OUT)

#### Plot During Simulation
plot.plots={'elateral':('sigma')}
plot.plot()

I also changed the recorder part. But when I run the code, I get the following Error:
File "<string>", line 1, in <module>
File "UCS.py", line 87, in recorder
NameError: global name 'elat_avg' is not defined

 AnnaMaria (annm5) said on 2019-11-28: #17

I am using more than one spaces
But when i copy paste my code in lanchpad it just removes them.
def recorder():
lateral()

 Jan Stránský (honzik) said on 2019-11-28: #18

def recorder():
elat_avg = lateral() # here