Metallic plate tension

Hi,

I did not get all your description (is this "plate" a sphere assembly ? that's it ? Is there no elements (like rigid walls) bounding this assembly) nor all your questions (maybe there was too many ?), but here are few elements:

- TriaxialStressController should work as good in tension as in compression, provided you switch the sign of the parameters (tension should be > 0 in recent versions). However, this engine is designed for a "cell" delimited with 6 rigid walls.

- if your sphere sample is a monosized regular (cristallographic) one without non-spherical boundary elements , your 3rd proposal (applying force directly to the boundary spheres) should just work. If N is the number of boundary spheres (easy to define for such sample), just apply to each sphere an external force F / N.

Note that shortening and precising your question could only help to get other (better ?) answers. ;-)

Jerome

thanks for response Jerome

Yes, in the picture for better explanation i just draw the rigid model)) In Yade now i use monosized regular specimen without non-spherical boundary elements, in future model will consists spheres with different size.

But if i set external force as you said, i think the distervution of strain and stress in model not correct (red color - max values, blue - min values)

displacement - dihttp://s28.postimg.org/ooujic9rh/puc1.png

stress - http://s29.postimg.org/dw4qu3xw7/puc1_stress.png

strain - http://s13.postimg.org/6bbtsjfav/puc1_strain.png

The problem is that spheres on the boundary have the lesser values of strain and stress, but displacement for them is highest. Also i can send result for the TriaxialStressController

thanks for response Jerome

Yes, in the picture for better explanation i just draw the rigid model)) In Yade now i use monosized regular specimen without non-spherical boundary elements, in future model will consists spheres with different size.

But if i set external force as you said, i think the distervution of strain and stress in model not correct (red color - max values, blue - min values)

displacement - http://s28.postimg.org/ooujic9rh/puc1.png

stress - http://s29.postimg.org/dw4qu3xw7/puc1_stress.png

strain - http://s13.postimg.org/6bbtsjfav/puc1_strain.png

The problem is that spheres on the boundary have the lesser values of strain and stress, but displacement for them is highest. Also i can send result for the TriaxialStressController

oh i forgot, of course plate is a sphere assembly

oh i forgot, of course plate is a sphere assembly

OK.
However, as a side note, let me draw to your attention that within DEM local stresses, as probably TesselationWrapper strains (I am not familiar with it), depend on relative displacements between particles, rather than absolute displacements.

So, maybe it is premature to conclude there is something wrong with the results shown in your comment #3.

Hello Alexander,

I’m new in Yade

welcome :-)

whenever you have a problem, send also a code reproducing your problem
(thanks for sending it without prior notification :-). However, try to send
it like MWE (minimal working example). It should be:
- minimal = short (it is difficult to find mistakes in long pieces of code
and it is usually much more annoying).
- working = before sending, test it to be sure it is working (it has no
syntax errors etc.).
together with above information, send also version of Yade you are using. I
tried your script, but it ended with segmentation fault. It could be
because I am using different version than you and if I knew your version, I
could try it using the correct one.

E.g. in the example you posted (next time) you could choose just one 'mode'
and shorten the saving (see below)

Also in one question, try to concentrate on one thing. Maybe (again next
time) I would ask one question on "what approach is the best for this
simulation" and another "I used this approach, why are the results so
strange?"

So it was some introduction into asking :-)

Now to your problem:
- when showing pictures with stress or strain, please be more specific on
what is plotted (some specific component? some invariant? something else?)
- note that 'young' and 'poisson' are not directly related to macroscopic
Young's modulus and Poisson's ratio. You can find some info in other
questions or documentation.
- could you please specify, what boundary conditions were used for the
pictures you sent?

One of the problem sources could be, that the simulation is run too shortly
and the state is far from quasi-static equilibrium and what you see is some
effect of inertial forces. It can be controlled using unbalancedForce()
function [4], something like (and deoending on the type of your simulation
control):

while unbalancedForce()>1e-2:
  O.run(100,True)

Some notes to the code (not crucial now):

when using interaction radius, it should be reset after the first step [3]

spheres=[]
> for i in range(0, 16):
> for j in range(0, 16):
> for k in range(0, 2):
> # for spheres which position satisfies condition below set constraint
> material
> if (i == 4 or i == 5 or i == 10 or i == 11) and (j == 4 or j == 5 or j ==
> 10 or j == 11):
> id =
> O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat_mod',radius=r))
> spheres.append(O.bodies[id])
> else:
> id =
> O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat_mod',radius=r))
> spheres.append(O.bodies[id])

you can use predefined function for this [1]:

from yade import pack
spheres = pack.regularOrtho(pack.inAlignedBox((0,0,0),(16.1,16.1,2.1)),r,0)
O.bodies.append(spheres)

I didn't use VTKExporter() because at a later stage I'll need another file
> extension to save the result data.

Just FYI (even if you will not use it in the future), to save spheres into
vtk file, you can also use predefined function [2]:

from yade import export
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[('velocity','b.state.vel'),('displacement','b.state.displ()'),('young','b.mat.young')])
# similarly for stress and strain

good luck using Yade :-)
cheers
Jan

[1] https://yade-dem.org/doc/yade.pack.html#yade.pack.regularOrtho
[2] https://yade-dem.org/doc/yade.export.html#yade.export.VTKExporter
[3] https://yade-dem.org/doc/user.html#creating-interactions
[4] https://yade-dem.org/doc/yade.utils.html#yade._utils.unbalancedForce

2015-07-09 21:26 GMT+02:00 Jérôme Duriez <
<email address hidden>>:

> Question #269063 on Yade changed:
> https://answers.launchpad.net/yade/+question/269063
>
> Status: Open => Answered
>
> Jérôme Duriez proposed the following answer:
> OK.
> However, as a side note, let me draw to your attention that within DEM
> local stresses, as probably TesselationWrapper strains (I am not familiar
> with it), depend on relative displacements between particles, rather than
> absolute displacements.
>
> So, maybe it is premature to conclude there is something wrong with the
> results shown in your comment #3.
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi, Jerome
Thank’s a lot for such a detail answer! So I’ll try to be more precise))
 So in the first picture (without spheres) I just draw the orthogonal plate like a real solid body to show what type of problem I have (without any connection to solving it with DEM).
As I said before boundary conditions are following
1. On face ABCD applied tensile force F parallel to OY axis (red arrow). it is necessary to say that force F is set for the whole face and not to specific points.
2. Face KGHO is fixed (Red dots)
So it means that the plate will stretch after force F has been applied to the side ABCD, because of the opposite side KGHO is fixed.
Objective is to compute distribution of strain, stress and displacement values on the plate after force has been applied. The final goal of my research is to solve this problem with 2 methods DEM and FEM and compare results.
In DEM I image it to do like this:
1) Represent plate like a set of spheres
2) Find correct way to set boundary conditions and apply them on spheres
3) Find strain, stress and displacement values for each sphere which represent the plate and save result data to file. (So now I save each sphere and it’s strain, stress, displacement attributes to VTK file which can be analyzed in Paraview)
Below I posted MWE example as you name it)) where permanent force is applied for spheres connected to boundary ABCD, spheres connected to KGHO is set to be fixed. (Thank’s for advising regularOrtho() function but still use simple loops:)
Yade version: 2015-03-17.git-16ecad0
////
###################################################
# define materials and model configuration

E = 1161e9 # Young's modulus of model
v = 0.33 # Poisson's ratio
p = 150e6 # initial stress value (positive - tension, negative - compression)

e = 0.02 # loading rate (strain rate) (positive - tension, negative - compression)
r = 0.5 # spheres radius

# Enlarge interaction radius between spheres using "interaction_radius" parameter (for example in uniax.py this value is 1.5)
interaction_radius = 1.5

# define plate material, create "dense" packing by setting friction to zero initially
O.materials.append(CpmMat(young=E,
           frictionAngle=0,
        poisson=v,
        density=4800,
        sigmaT=3.5e6,
        epsCrackOnset=1e-4,
        isoPrestress=0,
        relDuctility=30,
        label = 'mat_mod'))

# represent plate like a set of regular monosized set of spheres
# also set boundary conditions via predefined tensile force for spheres on ABCD and
# fixed spheres on KGHO
spheres=[]
for i in range(0, 16):
   for j in range(0, 16):
      for k in range(0, 2):
        id = O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat_mod',radius=r))
        spheres.append(O.bodies[id])
        if j == 15:
           O.forces.addF(id,(0,p,0),permanent=True) # Add force for all spheres connected to ABCD
        if j == 0:
           spheres[id].state.blockedDOFs='y' # Fixed all spheres connected to KGHO

###################################################
# define engines

# simulation loop
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius)]),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius)],
    [Ip2_CpmMat_CpmMat_CpmPhys()],
    [Law2_ScGeom_CpmPhys_Cpm()]
 ),
 CpmStateUpdater(realPeriod=1),
 NewtonIntegrator(damping=0.4)
]

###################################################
# start simulation and compute strain and stress

# try to run script with qt graphical interface
try:
   yade.qt.Controller(), yade.qt.View()
except:
   print 'Qt graphical interface is not avaliable'

# set the integration timestep to be 1/2 of the "critical" timestep
O.dt=.5*utils.PWaveTimeStep()

# compute strain via tesselation wrapper.
TW=TesselationWrapper()

# store current positions before simulation
TW.setState(0)

# run simulation
O.run(100,1)

# store positions after simulation (deformed state)
TW.setState(1)

# compute deformation for each body
TW.computeDeformations()

# compute stress tensor for each body
stresses = bodyStressTensors()

###################################################
# save data to vtk. file

n = len(spheres) # get number of particles
f = open('result.vtk','w')

# header
f.write('# vtk DataFile Version 3.0.\n')
f.write('comment\n')
f.write('ASCII\n\n')

# center position
string = str(n)
f.write('DATASET POLYDATA\n')
f.write('POINTS ')
f.write(string )
f.write(' double\n')
for b in spheres:
   for i in range(0,3):
      value = b.state.pos[i]
      string = str(value)
      f.write(string)
      f.write(' ')
   f.write('\n')
f.write('\n')

# radius
string = str(n)
f.write('POINT_DATA ')
f.write(string)
f.write('\n')
f.write('SCALARS radius double 1\n')
f.write('LOOKUP_TABLE default\n')
for b in spheres:
   value = b.shape.radius
   string = str(value)
   f.write(string)
   f.write('\n')
f.write('\n')

# strain
f.write('TENSORS strain_tensor float\n')
for b in spheres:
   for i in range(0,3):
      for j in range(0,3):
         value = TW.deformation(b.id,i,j)
         string = str(value)
         f.write(string)
         f.write(' ')
      f.write('\n')
   f.write('\n')
f.write('\n')

# stress
f.write('TENSORS stress_tensor float\n')
for b in spheres:
   for i in range(0,3):
      for j in range(0,3):
         value = stresses[b.id][i,j]
         string = str(value)
         f.write(string)
         f.write(' ')
      f.write('\n')
   f.write('\n')
f.write('\n')

# displacement
f.write('VECTORS displacement double\n')
for b in spheres:
   for i in range(0,3):
      value = b.state.displ()[i]
      string = str(value)
      f.write(string)
      f.write(' ')
   f.write('\n')
f.write('\n')

////

This code should give the results shown in the pictures:
displacement - http://s28.postimg.org/ooujic9rh/puc1.png
stress - http://s29.postimg.org/dw4qu3xw7/puc1_stress.png
strain - http://s13.postimg.org/6bbtsjfav/puc1_strain.png

These images shows set of spheres representing the plate (variable “spheres” in code) where each sphere is drawn with color depending on value of one chosen attribute (strain, stress or disp)
(as I said above each sphere contains values of 3 attributes strain, stress and disp)

So I don’t like stress and strain distribution, because for example spheres connected to ABCD and KGHO have minimal stress values.

And you already said my questions yourself :)

"I used this approach, why are the results so strange?"

“what approach is the best for this simulation?”

Also I wanna say that it’s very important for me to solve this problem:)

Soon I will post a picture with FEM results to compare, because in our test we consider FEM like main method and verify DEM via FEM.

I’ll be very glad if you can help me, thank’s a lot again

with regards, Alexander

Oh i'm so sorry I meant Jan in comment 8 :)

Jerome thank you too, but Jan's answer is more constructive

And also Jan how can use VTKExporter() to save stress and strain, because of that values saved in the "stresses" and "TW.deformation()" arrays in my example code

And they are not spheres attributes in my code sample

And they are not spheres attributes in my code sample

hm i need to do FEM example first to see what result is there, because now there is nothing to compare:) but anyway Jan i'm waiting for your opinion and recommendations about all of this. Am i on the right way?

Hi Alexander,

thanks for more explanation. However, some parts are still not clear :-)

- concerning stress and strain, you can't plot just stress, you plot its
component or some other scalar representation? What do you plot in the
pictures? yy component?

if j == 0:
> spheres[id].state.blockedDOFs='y' # Fixed all spheres
> connected to KGHO

- I don't believe the results are from MWE you posted, since the
displacement of spheres belonging to KGHO face is not constant (the blue
does not have constant shade), should be 0. Could you please rerun the
simulation and/or add coordinate system and scale to the pictures?
Moreover, are all KGHO points fixed only in y direction? If yes, then you
get perfect uniaxial loading and you you don't need any FEM for comparison,
just one line paper-pencil equation :-)

- anyway it will be difficult to calibrate Poisson's ratio to be 0.33 [2].
Just as the reminder, young and poisson parameters are not the values of
macroscopic Young's modulus an Poisson's ratio.

- you still did not reset the aabbEnlargeFactor
and interactionDetectionFactor [1]

- since the system is discrete, stress and strain values are just some
approximation. Moreover, at boundaries (=all spheres in your example :-),
the values might be approximated worse and also the physics itself
(distribution of interactions) is different than in "bulk" spheres.

- to save stresses (you can construct analogously the strains):

import __builtin__ # module visible from all python
stresses = bodyStressTensors()
__builtin__.my_stresses = stresses # make stresses visible from all python.
Be careful not to overwrite some important __builtin__ variable
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[...,('stress','__builtin__.my_stresses[b.id]')]) #
just stresses would not be visible from export module,
__builtin__.my_stresses are visible

cheers
Jan

[1] https://yade-dem.org/doc/user.html#creating-interactions
[2]
https://yade-dem.org/w/images/6/64/Stransky2010-Macroscopic-elastic-properties-of-particle-models.pdf

2015-07-10 21:21 GMT+02:00 Alexander <email address hidden>:

> Question #269063 on Yade changed:
> https://answers.launchpad.net/yade/+question/269063
>
> Alexander posted a new comment:
> hm i need to do FEM example first to see what result is there, because
> now there is nothing to compare:) but anyway Jan i'm waiting for your
> opinion and recommendations about all of this. Am i on the right way?
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hello Jan, I'm glad to see your answer.

Ok let's see consider the stress first.

At the previous picture i showed "magnitude" of stress (I can select such a mode in Paraview, but frankly i don't understand yet how it represents 3x3 matrix tensor http://i11.pixs.ru/storage/5/9/2/pic1JPG_5715109_18008592.jpg :)

Hm, I think the displacement of spheres belonging to KGHO is constant because of the minimal value is set to 0. (http://i11.pixs.ru/storage/6/0/7/pic2JPG_9158907_18008607.jpg). You are right, spheres of KGHO fixed just along Y direction. Also it seems to be a simple uniaxial loading, but as i mentioned in future i will replace specimen with another plate which consists inclusions of different metal (mo). Now i need to calibrate test that's why using simple uniform plate with stiffer material.

Today i created simple tension test in Ansys with the steel plate, the following picture shows stress distribution. (http://i11.pixs.ru/storage/6/3/3/pic3JPG_3955405_18008633.jpg). I want to see something like this in Yade, but after Young's modulus and Poisson's ratio were changed, the results become even worst. The spheres under constant force simply rip off the model. (http://i11.pixs.ru/storage/6/4/2/pic5JPG_4026499_18008642.jpg , huh i found how to show coordinate axis)

With you advise i'm proccessing to cut my MWE example. Now i try to use VTKExporter to save data but _builtin__.my_stresses doesn't work. That's why i save just dicplacement, you can uncomment variant with __builtin__.my_stresses . Pleace could you test it and said what's wrong with _builtin__.my_stresses.

/////

import __builtin__
from yade import export

###################################################
# define materials and model configuration

E = 2e11 # Young's modulus of model
v = 0.3 # Poisson's ratio
p = 150e6 # initial stress value (positive - tension, negative - compression)
d = 7850 # density

e = 0.02 # loading rate (strain rate) (positive - tension, negative - compression)
r = 0.5 # spheres radius

# Enlarge interaction radius between spheres using "interaction_radius" parameter (for example in uniax.py this value is 1.5)
interaction_radius = 1.5

# define plate material, create "dense" packing by setting friction to zero initially
O.materials.append(CpmMat(young=E,
           frictionAngle=0,
        poisson=v,
        density=d,
        sigmaT=3.5e6,
        epsCrackOnset=1e-4,
        isoPrestress=0,
        relDuctility=30,
        label = 'mat'))

# represent plate like a set of regular monosized set of spheres
# also set boundary conditions via predefined tensile force for spheres on ABCD and
# fixed spheres on KGHO
spheres=[]
for i in range(0, 16):
   for j in range(0, 16):
      for k in range(0, 2):
        id = O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat',radius=r))
        spheres.append(O.bodies[id])
        if j == 15:
           O.forces.addF(id,(0,p,0),permanent=True) # Add force for all spheres connected to ABCD
        if j == 0:
           spheres[id].state.blockedDOFs='y' # Fixed all spheres connected to KGHO

###################################################
# define engines

# simulation loop
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius,label='bo1s')]),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius,label='ig2ss')],
    [Ip2_CpmMat_CpmMat_CpmPhys()],
    [Law2_ScGeom_CpmPhys_Cpm()]
 ),
 CpmStateUpdater(realPeriod=1),
 NewtonIntegrator(damping=0.4)
]

###################################################
# start simulation and compute strain and stress

# try to run script with qt graphical interface
try:
   yade.qt.Controller(), yade.qt.View()
except:
   print 'Qt graphical interface is not avaliable'

# set the integration timestep to be 1/2 of the "critical" timestep
O.dt=.5*utils.PWaveTimeStep()

# compute strain via tesselation wrapper.
TW=TesselationWrapper()

# store current positions before simulation
TW.setState(0)

# run one single step
O.step()

# reset interaction radius to the default value
bo1s.aabbEnlargeFactor=1.0
ig2ss.interactionDetectionFactor=1.0

# run simulation
O.run(500,1)

# store positions after simulation (deformed state)
TW.setState(1)

# compute deformation for each body
TW.computeDeformations()

# compute stress tensor for each body
stresses = bodyStressTensors()

###################################################
# save data to vtk. file

__builtin__.my_stresses = stresses
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()')])
# with the command below script crashes!!!!
#vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','float (__builtin__.my_stresses[b.id])')])

////

So and the main question still the same, how can i simulate something like in Ansys picture with Yade. You also can test file result.vtk in Paraview with Glyph filter. I hope u will help me because i see u are professional in this area :). I will read u article but for me now to make working example is also very important.

with regards, Alexander

Hm, u right in picture http://i11.pixs.ru/storage/6/0/7/pic2JPG_9158907_18008607.jpg displacement of spheres belonging to KGHO it's not real 0, but may be it's because of rotation

Also how to represent data in paraview is here https://yade-dem.org/doc/user.html?highlight=paraview

Hi Alexander,

> At the previous picture i showed "magnitude" of stress (I can select such
> a mode in Paraview, but frankly i don't understand yet how it represents
> 3x3 matrix tensor
> http://i11.pixs.ru/storage/5/9/2/pic1JPG_5715109_18008592.jpg :)
>

exactly, it is always good to know what you are plotting :-) I don't know
what magnitude for tensors is, therefore I don't use it :-) paraview
numbering is 0=xx, 1=xy, 2=xz, 3=yx, 4=yy, 5=yz, 6=zx, 7=zy, 8=zz. So for
the first idea, I would use component 0 (xx)

>
> Hm, I think the displacement of spheres belonging to KGHO is constant
> because of the minimal value is set to 0.
> (http://i11.pixs.ru/storage/6/0/7/pic2JPG_9158907_18008607.jpg).

Again, here you plot magnitude, i.e. length of the vector. I would use 0
component (x direction)

> You are
> right, spheres of KGHO fixed just along Y direction. Also it seems to
> be a simple uniaxial loading, but as i mentioned in future i will
> replace specimen with another plate which consists inclusions of
> different metal (mo). Now i need to calibrate test that's why using
> simple uniform plate with stiffer material.
>

ok. What is actually your motivation to use DEM for this example?

>
> Today i created simple tension test in Ansys with the steel plate, the
> following picture shows stress distribution.
> (http://i11.pixs.ru/storage/6/3/3/pic3JPG_3955405_18008633.jpg). I want
> to see something like this in Yade, but after Young's modulus and
> Poisson's ratio were changed, the results become even worst. The spheres
> under constant force simply rip off the model.
>

Probably it is because you use CpmMat, intended for simulating concrete
including damage. To remain in elastic range, add neverDamage=True to
CpmMat definition.

> (http://i11.pixs.ru/storage/6/4/2/pic5JPG_4026499_18008642.jpg , huh i
> found how to show coordinate axis)
>
> With you advise i'm proccessing to cut my MWE example. Now i try to use
> VTKExporter to save data but _builtin__.my_stresses doesn't work. That's
> why i save just dicplacement, you can uncomment variant with
> __builtin__.my_stresses . Pleace could you test it and said what's wrong
> with _builtin__.my_stresses.
>

each component of stresses is tensor (Matrix3), you can't pass it to
float() function (Yade should give you TypeError). Just delete float and it
should work.

Hm, u right in picture
> http://i11.pixs.ru/storage/6/0/7/pic2JPG_9158907_18008607.jpg
> displacement of spheres belonging to KGHO it's not real 0, but may be
> it's because of rotation

it is because you plot magnitude, taking into account also lateral
displacements (which are not blocked and therefore nonzero).

cheers
Jan

Hi, Jan.

I'm sorry

vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','float (__builtin__.my_stresses[b.id])')]) is already my experiment,

because after

vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','__builtin__.my_stresses[b.id]')])

i got an error too.

About displacements u are also right, it's because yz - components are not fixed and i use magninute.

The rest of the questions i will answer later

with regards, Alexander

Hi Alexander,

>
> vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','float
> (__builtin__.my_stresses[b.id])')]) is already my experiment,
>
> because after
>
>
> vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','__builtin__.my_stresses[
> b.id]')])
>
> i got an error too.
>
>
If you get an error, please always put it to the email. It will help us to
understand more what is going on.

What version of Yade do you use?

cheers
Jan

Hi Jan,

My Yade version is

Yade version: 2015-03-17.git-16ecad0

What do u mean about e-mail: what adress of e-mail and what should i send? (code sample, error mesage..)

Also what should i use instead of CpmMat for such simulation?

With our colleagues i want to show that DEM is also applicable for such type of simulation, because it' much more simple to create calculational model in DEM. I told u at final step plate model will be more sophisticated, with inclusions of other more rigid material.

with regards, Alexander

Hi Alexander,

> What do u mean about e-mail: what adress of e-mail and what should i
> send? (code sample, error mesage..)
>

I meant that if you ask a question involving an error, put the error
message in the question :-) like

My error is:
Traceback (most recent call last):
  File "/home/honzik/bin/yade-trunk", line 182, in runScript
    execfile(script,globals())
  File "a.py", line 104, in <module>

vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','float
(__builtin__.my_stresses[b.id])')])
  File "/home/honzik/lib/x86_64-linux-gnu/yade-trunk/py/yade/export.py",
line 431, in exportSpheres
    test = eval(command) # ... eval one example to see what type (float,
Vector3, Matrix3) the result is ...
  File "<string>", line 1, in <module>
NameError: name '__builtin__' is not defined

>
> Also what should i use instead of CpmMat for such simulation?
>

it depends on what you want to simulate. If only elastic behavior, then
CpmMat (with neverDamage=True) is fine. You can also try CohFrictMat.

>
> With our colleagues i want to show that DEM is also applicable for such
> type of simulation,

definitely it is.

> because it' much more simple to create calculational
> model in DEM.

I would be careful with terms like "much more simple" :-)

concerning __builtin__ issue, sorry for my mistake, the idea is this:

import __builtin__
...
stresses = bodyStressTensors
__builtin__.my_stresses = stresses
...
vtk.exportSpheres(what=[...,('stress','my_stresses[b.id]')]) # my_stresses
*without* __builtin__.

Once more, be careful using __builtin__ module not to overwrite some
important variable. E.g.

print tuple([1,2,3]) # no problem
import __builtin__
__builtin__.tuple = dict
print tuple([1,2,3]) # some strange errors

cheers
Jan

Hi, Jan.

thank's for fixing __builtin__.my_stresses. Now i'm wondering how to import strains because of TW.deformation(b.id,i,j) works only for one sphere, i think i need to create new array strains[] and put all values in it first and then use __builtin__, or it can be done more elegantly?

Now the problem itself:

1) I'v tried to fix all u remarks (MWE is in the end of the message).

2) Well, let's consider the first component of displacement which i get with Yade - (http://i11.pixs.ru/storage/5/5/3/pic1JPG_1997634_18020553.jpg) and Ansys - (http://i11.pixs.ru/storage/5/5/7/pic3JPG_1380016_18020557.jpg)

As u can see they are different, Yade has values' range [-0.1038,0.1038] Ansys range is [-0,0018636,0,0018636]

3) The same thing with the x component of normal stress (first component of stress tensor)
Yade - (http://i11.pixs.ru/storage/5/5/9/pic2JPG_1013994_18020559.jpg) Ansys - (http://i11.pixs.ru/storage/5/6/4/pic4JPG_1715017_18020564.jpg) Yade range = [-3.32165e7,6.07687e7] Ansys range = [-2,2244e7,9,8264e7].

The size of model in Ansys is the same (16,16,2) - (http://i11.pixs.ru/storage/5/8/8/pic6JPG_8207000_18020588.jpg) also u can see computational mesh (http://i10.pixs.ru/storage/5/9/1/pic7JPG_2524627_18020591.jpg)

So, any ideas?

MWE

///

import __builtin__
from yade import export

###################################################
# define materials and model configuration

E = 2e9 #11 # Young's modulus of model
v = 0.3 # Poisson's ratio
p = 150e6 # initial stress value (positive - tension, negative - compression)
d = 7850 # density

e = 0.02 # loading rate (strain rate) (positive - tension, negative - compression)
r = 0.5 # spheres radius

# Enlarge interaction radius between spheres using "interaction_radius" parameter (for example in uniax.py this value is 1.5)
interaction_radius = 1.5

# define plate material, create "dense" packing by setting friction to zero initially
O.materials.append(CpmMat(young=E,
           frictionAngle=0,
        poisson=v,
        density=d,
        sigmaT=3.5e6,
        epsCrackOnset=1e-4,
        neverDamage=True,
        isoPrestress=0,
        relDuctility=30,
        label = 'mat'))

# represent plate like a set of regular monosized set of spheres
# also set boundary conditions via predefined tensile force for spheres on ABCD and
# fixed spheres on KGHO
spheres=[]
for i in range(0, 16):
   for j in range(0, 16):
      for k in range(0, 2):
        id = O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat',radius=r))
        spheres.append(O.bodies[id])
        if j == 15:
           O.forces.addF(id,(0,p,0),permanent=True) # Add force for all spheres connected to ABCD
        if j == 0:
           spheres[id].state.blockedDOFs='xyzXYZ' # Fixed all spheres connected to KGHO

###################################################
# define engines

# simulation loop
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius,label='bo1s')]),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius,label='ig2ss')],
    [Ip2_CpmMat_CpmMat_CpmPhys()],
    [Law2_ScGeom_CpmPhys_Cpm()]
 ),
 CpmStateUpdater(realPeriod=1),
 NewtonIntegrator(damping=0.4)
]

###################################################
# start simulation and compute strain and stress

# try to run script with qt graphical interface
try:
   yade.qt.Controller(), yade.qt.View()
except:
   print 'Qt graphical interface is not avaliable'

# set the integration timestep to be 1/2 of the "critical" timestep
O.dt=.5*utils.PWaveTimeStep()

# compute strain via tesselation wrapper.
TW=TesselationWrapper()

# store current positions before simulation
TW.setState(0)

# run one single step
O.step()

# reset interaction radius to the default value
bo1s.aabbEnlargeFactor=1.0
ig2ss.interactionDetectionFactor=1.0

# run simulation, until static equilibrium will not reached
while unbalancedForce()>1e-2:
   O.run(10,True)

# store positions after simulation (deformed state)
TW.setState(1)

# compute deformation for each body
TW.computeDeformations()

# compute stress tensor for each body
stresses = bodyStressTensors()

###################################################
# save data to vtk. file
__builtin__.my_stresses = stresses
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','my_stresses[b.id]')])

///

with regards, Alexander

Pressure F on ABCD is also the same: http://i11.pixs.ru/storage/6/5/1/pic8JPG_6103973_18020651.jpg

Hi Alexander,

> thank's for fixing __builtin__.my_stresses. Now i'm wondering how to
> import strains because of TW.deformation(b.id,i,j) works only for one
> sphere, i think i need to create new array strains[] and put all values
> in it first and then use __builtin__, or it can be done more elegantly?
>
>
I don't know TessalationWrapper, but constructing the values can be done
relatively elegantly:

strains = [Matrix3(*[TW.deformation(b.id),i,j) for i in xrange(3) for j in
xrange(3)]) for b in O.bodies]

>
> 2) Well, let's consider the first component of displacement which i get
> with Yade -
> (http://i11.pixs.ru/storage/5/5/3/pic1JPG_1997634_18020553.jpg) and
> Ansys - (http://i11.pixs.ru/storage/5/5/7/pic3JPG_1380016_18020557.jpg)
>
> As u can see they are different, Yade has values' range
> [-0.1038,0.1038] Ansys range is [-0,0018636,0,0018636]
>
> 3) The same thing with the x component of normal stress (first component
> of stress tensor)
> Yade - (http://i11.pixs.ru/storage/5/5/9/pic2JPG_1013994_18020559.jpg)
> Ansys - (http://i11.pixs.ru/storage/5/6/4/pic4JPG_1715017_18020564.jpg)
> Yade range = [-3.32165e7,6.07687e7] Ansys range = [-2,2244e7,9,8264e7].
>
> The size of model in Ansys is the same (16,16,2) -
> (http://i11.pixs.ru/storage/5/8/8/pic6JPG_8207000_18020588.jpg) also u
> can see computational mesh
> (http://i10.pixs.ru/storage/5/9/1/pic7JPG_2524627_18020591.jpg)
>
> So, any ideas?
>

The difference is most likely because 'young' and 'poisson' arguments of
CpmMat are interaction parameters, *not* the macroscopic paremters (e.g.
different material in FEM and DEM). A standard way is to calibrate these
micro parameters to get macroscopically what you want. In one of my
previous messages you can find a reference on this.

The perfect uniaxial condition (recall the paper-pencil example I suggested
:-) may be a good starting point. Then you would have constant stress and
stress, linear displacement, Young's modulus and Poisson's ratio,
everything computable by hand. Also for the first FEM and DEM comparison, I
would use uniaxial tension (most simple loading scenario) instead of that
you are using now..

Concerning the plotting, I messed a bit the axes, what I meant was to print
both displacement and stresses in the direction of applied force. Anyway,
if you want the same results, you have to play with input material
parameters :-)

cheers
Jan

Hi, Jan.

First thanks for strains = [Matrix3(*[TW.deformation(b.id),i,j) for i in xrange(3) for j in
xrange(3)]) for b in O.bodies] i will check it.

Now the problem, so i also messed with value of Young's module taken 2e9 instead of 2e11 :), that's why the results was so different, now they are a little bit closer , please see pictures along the direction of pressure F.

displacement Yade - http://i11.pixs.ru/storage/5/9/0/pic1JPG_8725641_18032590.jpg Ansys - http://i11.pixs.ru/storage/6/1/0/pic4JPG_2632224_18032610.jpg

stress Yade - http://i11.pixs.ru/storage/7/6/4/pic2JPG_7240037_18032764.jpg Ansys - http://i11.pixs.ru/storage/6/5/6/pic5JPG_6182930_18032656.jpg.

Talking about stress besides values range i also don't like distribution because of the spheres on ABCD and KGOH has minimal values (marked with red rectangules) . Stress distribution should be something like in Ansys, because we checked it it with another FEM simulator (http://i11.pixs.ru/storage/7/5/5/pic7JPG_3148063_18032755.jpg).

You already made comment: "at boundaries (=all spheres in your example :-), the values might be approximated worse and also the physics itself (distribution of interactions) is different than in "bulk" spheres." So should i left this situation like that or do something? Anyway showing stress distribution like that is not ok.

Well now about calibrating 'young' and 'poisson' arguments.

 I fluently read your article (soon I'll try read it in details). As i understand the final relation between macroscopic material properties and microscopic model parameters described by the formulas (highlighted in red)(http://i11.pixs.ru/storage/8/1/0/pic3JPG_6512680_18032810.jpg) .

So should i use these formulas or set 'young' and 'poisson' empirically (by changing their values and compare with Ansys)?

If i need to use formulas the questions are:

1) where can i get k_N and K_S params in CpmMat?
2) Also how to get L_c, as i understand L_c is the distance between centers of each 2 sphere bounded by link (so called bars length),
3) Is it applicable only for model consists of rigid spheres with uniform radius R?

If i will try second variant (i think it's simplier :) what is the best stategy to do it? Because i have 2 parameters and should change it separately (may be i should fix one and change another, or for example change both of them at one step)?

with regards, Alexander

I want to add about "stress distribution". May be i should use larger quantity of spheres in packing or do something with pressure F applied on ABCD spheres at the end..ets?

Hi, Jan

As you recommend I computed deformation along OY axis by hand; it's around FL/ES=PSL/ES=PL/E=(150e6 )*16/2e11= 150*16/2e5= (15*16)/(2*10000)=0,012 where P is pressure like in my example. So the result is like in Ansys http://i11.pixs.ru/storage/6/5/6/pic5JPG_6182930_18032656.jpg.

with regards, Alexander

Oh meant this http://i11.pixs.ru/storage/6/1/0/pic4JPG_2632224_18032610.jpg

So and if we set simple uniaxial tension in Ansys the results became exact

disp - http://i11.pixs.ru/storage/4/7/5/pic1JPG_5894468_18044475.jpg

strain = 0,012 / 16 = 0,00075 - http://i11.pixs.ru/storage/4/8/6/pic2JPG_4441029_18044486.jpg

stress = http://i11.pixs.ru/storage/5/0/2/pic3JPG_5210891_18044502.jpg

By modeling uniaxial tension with Yade I achieve same displacement with young = E = 1.338e11, v = 0.3. (http://i11.pixs.ru/storage/6/2/1/pic4JPG_8549414_18044621.jpg) but the values range and distribution of stress are not correct at all (http://i11.pixs.ru/storage/6/3/8/pic5JPG_4880801_18044638.jpg).

So any ideas Jan what i gotta do?

with regards, Alexander

By modeling uniaxial tension with Yade I achieve same displacement with young = E = 1.338e11, v = 0.3. (http://i11.pixs.ru/storage/6/2/1/pic4JPG_8549414_18044621.jpg) but the values range and distribution of stress are not correct at all (http://i11.pixs.ru/storage/6/3/8/pic5JPG_4880801_18044638.jpg).

So any ideas Jan what i gotta do?

with regards, Alexander

Oh forget to choose YY normal component http://i11.pixs.ru/storage/6/7/5/pic5JPG_8836063_18044675.jpg but result is the same

Also code reproduces last picture

import __builtin__
from yade import export

###################################################
# define materials and model configuration

E = 1.338e11#2e11 # Young's modulus of model
v = 0.3#0.3 # Poisson's ratio
p = 150e6 # initial stress value (positive - tension, negative - compression)
d = 7850 # density

r = 0.5 # spheres radius

# Enlarge interaction radius between spheres using "interaction_radius" parameter (for example in uniax.py this value is 1.5)
interaction_radius = 1.5

# define plate material, create "dense" packing by setting friction to zero initially
O.materials.append(CpmMat(young=E,
           frictionAngle=0,
        poisson=v,
        density=d,
        sigmaT=3.5e6,
        epsCrackOnset=1e-4,
        neverDamage=True,
        isoPrestress=0,
        relDuctility=30,
        label = 'mat'))

# represent plate like a set of regular monosized set of spheres
# also set boundary conditions via predefined tensile force for spheres on ABCD and
# fixed spheres on KGHO
spheres=[]
for i in range(0, 16):
   for j in range(0, 16):
      for k in range(0, 2):
        id = O.bodies.append(sphere([i+0.5,j+0.5,k+0.5],material='mat',radius=r))
        spheres.append(O.bodies[id])
        if j == 15:
           O.forces.addF(id,(0,p,0),permanent=True) # Add force for all spheres connected to ABCD
        if j == 0:
           #spheres[id].state.blockedDOFs='xyzXYZ' # Fixed all spheres connected to KGHO
           O.forces.addF(id,(0,-p,0),permanent=True)

###################################################
# define engines

# simulation loop
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius,label='bo1s')]),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius,label='ig2ss')],
    [Ip2_CpmMat_CpmMat_CpmPhys()],
    [Law2_ScGeom_CpmPhys_Cpm()]
 ),
 CpmStateUpdater(realPeriod=1),
 NewtonIntegrator(damping=0.4)
]

###################################################
# start simulation and compute strain and stress

# try to run script with qt graphical interface
try:
   yade.qt.Controller(), yade.qt.View()
except:
   print 'Qt graphical interface is not avaliable'

# set the integration timestep to be 1/2 of the "critical" timestep
O.dt=.5*utils.PWaveTimeStep()

# compute strain via tesselation wrapper.
TW=TesselationWrapper()

# store current positions before simulation
TW.setState(0)

# run one single step
O.step()

# reset interaction radius to the default value
bo1s.aabbEnlargeFactor=1.0
ig2ss.interactionDetectionFactor=1.0

# run simulation, until static equilibrium will not reached
while unbalancedForce()>1e-2:
   O.run(10,True)

# store positions after simulation (deformed state)
TW.setState(1)

# compute deformation for each body
TW.computeDeformations()

# compute stress tensor for each body
stresses = bodyStressTensors()

# get strain tensor for each body
strains = [Matrix3(*[TW.deformation(b.id,i,j) for i in range(0,3) for j in range(0,3)]) for b in O.bodies]

###################################################
# save data to vtk. file
__builtin__.my_stresses = stresses
__builtin__.my_strains = strains
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','my_stresses[b.id]'),('strains','my_strains[b.id]')])

with regards, Alexander

Short comments (I didn't follow the complete discussion, sorry):

1/ About comparison of stress:
- stress Yade - http://i11.pixs.ru/storage/7/6/4/pic2JPG_7240037_18032764.jpg Ansys - http://i11.pixs.ru/storage/6/5/6/pic5JPG_6182930_18032656.jpg.

If you compare the number of Gauss points in Ansys with the number of spheres in you DEM model I guess you'll find a large difference, hence the comparison is not fair. It suggests a more general question: why simulating an elastic plate with DEM at all?! With such regular array of spheres you get a very peculiar orthotropic behaviour which will never match the results obtained for isotropic materials.

2/ Micro-macro relation for Young's modulus in regular array is trivial: Emacro = Emicro (put a factor 0.5 maybe, not completely sure now). If your material type does not have Emicro there must be something equivalent.
Poisson's ratio is null in the directions of the regular array (for all others directions it is a different story...)

Hi Bruno.

Ok the mesh constructed in Ansys is also not very fine (http://i10.pixs.ru/storage/5/9/1/pic7JPG_2524627_18020591.jpg) i'll try more number spheres but i think it also wouldn't change the results around boundary (especially stress).

As i said before:

With our colleagues i want to show that DEM is also applicable for such type of simulation, because it' much more simple to create calculational model in DEM. At the final step plate model will be more sophisticated, with inclusions of other more rigid material.

About regularity: should i use random packing instead of regular array or what?

And the main question is still the same: How to get something similar to Ansys results with Yade? Let it be in the first step simple uniaxial tension of steel plate with the predefined pressure. Like in the last pictures. I mean deformation, stress and strains (i'm interesting in value range and distribution )

Yet Another Detailled Explanation (attempt, at least) about your "main question":

The mechanics behind the Ansys and Yade computations have *NOTHING* in common:

- discrete mechanics with (almost) spring constitutive laws (interaction forces <-> relative displacements) for Yade on one hand

- continuum mechanics with elastic or elasto-plastic constitutive relations (stress <-> strain) for Ansys on the other hand.

Obviously, since both methods are used sometimes to describe the same physical reality, it should be possible to let match the results of these two methods....
However, this requires a non-trivial calibration procedure, so that the parameters of the two models (Yade and Ansys) lead to similar results.
As a matter of fact, this is non-trivial at all, since Ansys parameters have nothing in common with Yade parameters: they do not enter in the same equations !

As for the outputs of Yade and Ansys models, you may again first consider that they have nothing in common ! Yade describes a particulate system (=> displacement, forces...), Ansys describes a continuum (=>stress and strain tensors).
Again, it is however possible to link to some extent these different quantities. But, after some reasoning !

For instance, stress tensor definitions -- from forces within a discrete system -- makes really sense only for high number of particles (=> Representative Elementary Volume notion). Which is not the case with getBodyStress(). You still might be obtain some insights from such function, but, in the general case, you will never obtain an uniform stress field simulating uniaxial tension as you may get with Ansys !

I understand that DEM and FEM are different methods, that why i still cannot solve my problem:) Ok I'll try getBodyStress() and more dense packing.

Also i don't expect perfect equality but anyway results gotta be more closer. I don't like rows of blue spheres along the boundary where pressure was apllied http://i11.pixs.ru/storage/6/3/8/pic5JPG_4880801_18044638.jpg. Also the range of stress values is very wide [1.48e8,3e8] in ansys it's 1,5e8.

May be i got to use another way to set pressure? Or some other class of material instead of CmpMat?

The best way will be give me some instructions in details what to do like:
1) Define material with the following class....(CmpMat e.g.)
2) To calibrate young and poisson values u should use following algorithm...
1) Create spheres packing (number of spheres, their relative location...etc)
2) Apply predefined pressure with ... function
3) Compute stresses like this...
etc..

Thank everybody for any help

with regards, Alexander

Hi Alexander,

> 1) Define material with the following class....(CmpMat e.g.)
>

for elastic behavior, any material (able to transfer tension) is fine. No
problem here

> 2) To calibrate young and poisson values u should use following
> algorithm...
>

For elastic material, you want to fit 2 macroscopic constants, e.g. Young's
modulus and Poissons's ratio. The main conclusion from the paper:
- The Young's modulus is directly proportional to 'young' CpmMat
parameters, so it can be fit trivially.
- Poisson's ratio does depend on 'poisson' CpmMat parameter, but
nonlinearly, so you can run simulations with different values to see the
relationship. As it is 1 variable, the optimization is not so difficult.

> 1) Create spheres packing (number of spheres, their relative
> location...etc)
>

Well, this is mainly up to you :-) Depending on your purposes (geometry,
loading), you can use regular packing, or random one
- regular packing is fine from my point of view for regular geometries.
- irregular is usually random so you will not get constant results. Also
you should use large amount of particles to get realistic results (In the
limit of infinite number of particles, you would get perfect results :-).
Moreover the boundaries will not be smooth anymore
- you may consider using some different discrete approach, something like
Cusatis' lattice model. You have very regular interaction network, so you
can also solve one cell of the truss-like interaction network analytically
"by hand".

> 2) Apply predefined pressure with ... function
>

this way (prescribing force on each particle from boundary layer) is the
most realistic

> 3) Compute stresses like this...
>

probably the biggest problem is here :-) If your displacements corresponds
to macroscopic Young's modulus and Poisson's ratio, then they are correct.
The stress and strain is just some post-processing and approximation from
discrete forces and displacements, so perfect results are not expectable :-)

Two points here:
- apart from stress, plot also contact forces. You will see, that at the
boundary (parallel to force) layer, the forces are different. It is because
the interactions distribution is different than "in the middle" (same
interactions in parallel direction, but only half in perpendicular and
diagonal direction). So if you are happy with middle stresses and strains,
you are on the right way :-)
- you can also put lower force value for conrner particles, as the force is
transmitted by this nonstandard parallel boundary layer. Do not forget to
redistribute this jump to the "middle" forces to preserve total force value
:-) You can again try different "jump" values to see the response, and
choose the correct "jump" value. The same approach you would choose, if you
directly apply nodal forces on FEM mesh and you want to get constant
stress/strain..
- Stress is computed from interactions. The fixed boundary and boundary
with prescribed force has half parallel interactions, so (very rouhly) you
get half stress. See [1], you can also compute it analytically by hand :-)
- The more particles you use, the smaller would be this boundary effects on
the overall results
- finally, you can try lower value of unbalanced force

cheers
Jan

PS: @Bruno: there are also diagonal bonds, so the interaction network is
relatively dense and the isotropy is moreless ok

[1] https://yade-dem.org/doc/yade.utils.html#yade._utils.bodyStressTensors

>the mesh constructed in Ansys is also not very fine

The real question is the number of nodes/gauss points in each element times the number of elements, not the number of elements alone. Looking at you stress field, I guess they are high order elements (i.e. many DOFs per element).

>About regularity: should i use random packing instead of regular array or what?

Yes certainly if you want to simulate an isotropic material. The problem if you do that is that per-body stress tensors will become highly heterogeneous (while it is ok to use it in the regular packing you have now), then you will need spatial averaging of some sort and sufficiently high discretization level for those local averages to be doable in small regions (compared to problem size).

>And the main question is still the same: How to get something similar to Ansys results with Yade?

To me it is still not very clear what/why you think is wrong after a quick overview of all above messages. Basically, the solution (displacement and forces - converted into strain/stress if you like) for such regular array can be derived with paper and pen, so I would suggest to divide in two sub-problems:
P1- Does the analytical match ANSYS (certainly not since a sphere array is not isotropic)?
P2- Does DEM match the analytical expressions? If not there can be a mistake in the definition of boundary condition, a non-static state, or something like that.
For the moment it is not clear if you struggle with P1, P2, or both of them. It makes the problem solving difficult.
I understand you'd like a step by step answer on how to do things but there is no well posed problem yet in my opinion.

After reading the last messages I recognize I'm still a few steps behind Jan, so you can waste my reply maybe. :)

@Jan
Ok for diagonal interactions, I didn't notice. It is indeed more isotropic, but on another hand isn't it what explains the special trend of stress on free boundaries? This boundary effect is still a problem in my eyes (half stress where force is applied is ok, it is a post processing artifact).

Bruno

Last thing, as far as per-body stress is a topic of concern: getBodyStress is exact.
It should be exactly equal to applied stress as long as you account for the volume occupied each sphere consistently (in a porous medium the local solid phase stress is larger than the macroscale stress).
In your regular array, you can apply a simple conversion rule for each sphere:
macroStress = sphereStress*sphereVolume/cubeVolume
where the cube circumbscribes the sphere.
The average of the above (in a row or in the full domain) must be exactly equal to the ansys stress - which is itself exactly equal to applied stress.

Hello Jan and Bruno, thank you for your answers?

Now my problems:

1) By your advise I increase the dense of packing (still left it regular), but after this displacement of spheres is also changed.

Dislplacements along Y axis gotta be around [-0.006, 0.006]

Pay attention! In my simulation i applied force F like pressure.

Test 2 variants, applying force F directly to spheres and multiplying by sphere radius to take into account spheres size.
For each variant different value of young modulus was used and it's calibrated for coarse variant of packing

First variant: coarse - http://i11.pixs.ru/storage/7/5/8/pic1JPG_5947378_18056758.jpg fine - http://i11.pixs.ru/storage/7/6/8/pic2JPG_7134594_18056768.jpg

Second variant: coarse - http://i11.pixs.ru/storage/7/7/4/pic3JPG_5739457_18056774.jpg fine - http://i11.pixs.ru/storage/7/7/7/pic4JPG_2331364_18056777.jpg

You can see that desired displacement changed after both variants.

So is it correct to apply force (pressure) to sphere without taking into account it's size, or i need to multiply by it's radius, diameter or volume? And what should i do with changing displacement after different dense of packing was applied?

Well now let's discuss stress values distribution

In this picture i calibrated young to get appropriate displacement for the fine packing http://i11.pixs.ru/storage/8/3/6/pic5JPG_4896692_18056836.jpg

The values range of stress is here: http://i11.pixs.ru/storage/8/5/3/pic6JPG_4993977_18056853.jpg
If conversion macroStress = sphereStress*sphereVolume/cubeVolume will be used we get
something like 3,2321e8 / 2 = 1,69e8 which is similar to 1,5e8. (like exact variant in Ansys)

As you can see in the following pictures such maximum values have spheres inside the packing surrounded by one boundary layer of spheres.

whole model - http://i11.pixs.ru/storage/9/2/3/pic8JPG_2697509_18056923.jpg

slice along XOZ - http://i11.pixs.ru/storage/9/2/5/pic9JPG_7187525_18056925.jpg

slice along XOY http://i11.pixs.ru/storage/9/3/0/pic10JPG_8961910_18056930.jpg

So may be i should just create some additional boundary layer of spheres sorrounded the model and then delete it? Because your advise Jan:

>- you can also put lower force value for conrner particles, as the force is
>transmitted by this nonstandard parallel boundary layer. Do not forget to
>redistribute this jump to the "middle" forces to preserve total force value
>:-) You can again try different "jump" values to see the response, and
>choose the correct "jump" value. The same approach you would choose, if you
>directly apply nodal forces on FEM mesh and you want to get constant
>stress/strain..

couldn't help with spheres where force is applied, they have fixed applied force on them. (in the pictures they are always blue)

Here yours comment about that:
>- Stress is computed from interactions. The fixed boundary and boundary
>with prescribed force has half parallel interactions, so (very rouhly) you
>get half stress. See [1], you can also compute it analytically by hand :-)

Also Bruno about your formula macroStress = sphereStress*sphereVolume/cubeVolume. How it will work for non regular packing? Because in my final model there will be non regular packing areas inside the plate.

> Does the analytical match ANSYS (certainly not since a sphere array is not isotropic)?
Yes Ansys solution is exact for uniaxial tension

>P2- Does DEM match the analytical expressions? If not there can be a mistake in the definition of boundary condition, a non->static state, or something like that.
What do you mean here, the result with etalon stress and dicplacement in Ansys doesnt match. Honestly i'm still nor ok with DEM theory:)

Do you know some short and easily understood literature for it. For example at the beginning of Jan's article strain and the contact spheres by springs are described quite well, but there's now stress there. Also what type of contact used in Yade? I think the same like in Jan's article. May be you have something else Jan, especially for my type of simulation?

Also Jan i will read about https://yade-dem.org/doc/yade.utils.html#yade._utils.bodyStressTensors.

>Two points here:
>- apart from stress, plot also contact forces. You will see, that at the
>boundary (parallel to force) layer, the forces are different. It is because

How to plot contact forces and how much vectors of force will have each sphere?

> The real question is the number of nodes/gauss points in each element times the number of elements,
You are right that's why i tried more dense packing, but the results is not ok yet

>Ok for diagonal interactions, I didn't notice. It is indeed more isotropic, but on another hand isn't it what explains the special >trend of stress on free boundaries? This boundary effect is still a problem in my eyes (half stress where force is applied is ok, >it is a post processing artifact).
Here i think only Jan can help, he is much more professional than me:)

>For the moment it is not clear if you struggle with P1, P2, or both of them. It makes the problem solving difficult.
Of course i have struggle with P2, FEM is fine:) For uniaxial tension which i test now, Ansys has perfect results (you can see it in the previous comments on pictures)

>So if you are happy with middle stresses and strains,you are on the right way :-)
Jan, I'm not happy yet with middle stresses because they are bigger:) And by the way what do you think about Bruno's proposal to devide macroStress = sphereStress*sphereVolume/cubeVolume and what should i do in the non-regular packing?
As said in the future there will be non-regular areas inside.

I'm still waiting for help:)

with regards, Alexander

I'v already read Yade manual with DEM background but there is no simple formulas for strain and stress, for example in you article Jan strain described much better. Also it's interesting to read about interactions laws in CpmMat material model which i used now.

Also the last MWE:

import __builtin__
from yade import export

###################################################
# define materials and model configuration

E = 5.2e11#0.65e11#1.338e11#2e11 # Young's modulus of model need to be calibrated
v = 0.3# Poisson's ratio
P = 150e6 # initial pressure value
d = 7850 # density

r = 0.25 # spheres radius

# Enlarge interaction radius between spheres using "interaction_radius" parameter (for example in uniax.py this value is 1.5)
interaction_radius = 1.5

F = P*r

# define plate material, create "dense" packing by setting friction to zero initially
O.materials.append(CpmMat(young=E,
           frictionAngle=0,
        poisson=v,
        density=d,
        sigmaT=3.5e6,
        epsCrackOnset=1e-4,
        neverDamage=True,
        isoPrestress=0,
        relDuctility=30,
        label = 'mat'))

# represent plate like a set of regular monosized set of spheres
# also set boundary conditions via predefined tensile force for spheres on ABCD and
# fixed spheres on KGHO
O.bodies.append(pack.regularOrtho(pack.inAlignedBox((0,0,0),(16.1,16.1,2.1)),r,0))
for b in O.bodies:
    id = b.id
    if b.state.pos[1]+r > 15.9:
      O.forces.addF(id,(0,F,0),permanent=True)
    if b.state.pos[1]-r < 0.1:
      O.forces.addF(id,(0,-F,0),permanent=True)

###################################################
# define engines

# simulation loop
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius,label='bo1s')]),
 InteractionLoop(
 [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius,label='ig2ss')],
    [Ip2_CpmMat_CpmMat_CpmPhys()],
    [Law2_ScGeom_CpmPhys_Cpm()]
 ),
 CpmStateUpdater(realPeriod=1),
 NewtonIntegrator(damping=0.4)
]

###################################################
# start simulation and compute strain and stress

# try to run script with qt graphical interface
try:
   yade.qt.Controller(), yade.qt.View()
except:
   print 'Qt graphical interface is not avaliable'

# set the integration timestep to be 1/2 of the "critical" timestep
O.dt=.5*utils.PWaveTimeStep()

# run one single step
O.step()

# reset interaction radius to the default value
bo1s.aabbEnlargeFactor=1.0
ig2ss.interactionDetectionFactor=1.0

# run simulation, until static equilibrium will not reached
while unbalancedForce()>1e-2:
   O.run(10,True)

# compute stress tensor for each body
stresses = bodyStressTensors()

###################################################
# save data to vtk. file
__builtin__.my_stresses = stresses
vtk = export.VTKExporter('result')
vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','my_stresses[b.id]')])

Also i mean displacement in the terms of deformation of the model

Hello Alexander,

First of all, this question is becoming too complex, solving a lot of
topics at once. I would vote to close this one and open separate question
for each problem (stress computation, material model, ...)

1) By your advise I increase the dense of packing (still left it
> regular), but after this displacement of spheres is also changed.

sure, you get correct results in the limit case of infinte number of
particles. What make problem are the boundary particles, which have
different number of interactions than "middle" particles. In your first
attemt, ALL particles were the boundary particles :-) Now if you refine
even more, the difference should be smaller..

Pay attention! In my simulation i applied force F like pressure.

> Test 2 variants, applying force F directly to spheres and multiplying by
> sphere radius to take into account spheres size.
> For each variant different value of young modulus was used and it's
> calibrated for coarse variant of packing

 Here I do not understand at all. In DEM, there is nothing like pressure
(!!!!), only discrete forces. In your code you apply (0,F,0) force to each
particle:

> F = P*r
> ...

O.forces.addF(id,(0,F,0),permanent=True)

If you sum all forces on one face, you should get total force, like
pressure*area on physical sample. In case of regular packing, you should
use P*r*r, r*r being the aree represented by one particle (just by
comparing physical dimensions, multipliyng pressure by length does not
yield force)

couldn't help with spheres where force is applied, they have fixed
> applied force on them. (in the pictures they are always blue)

sure you could help, just use different formula or some correction :-)
bodyStressTensors have a few assumptions:
- the sphere is in static equilibrium (you probably fulfil this)
- there is no (non-interaction) force applied on the particle (which is not
the case of boundary particles with prescribed force)
- the particle have "enough" interactions, the more in all directions, the
better. Ok for "inner particle", not very ok for the boundary ones, as the
interactions are only in one halfspace

Also Bruno about your formula macroStress =
> sphereStress*sphereVolume/cubeVolume. How it will work for non regular
> packing? Because in my final model there will be non regular packing
> areas inside the plate.

The formula works for both regular and irregular packings, just have the
above assumptions in head. Also not another Bruno's comment that for
irregular case, the stress would be oscillating quantity

Do you know some short and easily understood literature for it. For example
> at the beginning of Jan's article strain and the contact spheres by springs
> are described quite well, but there's now stress there. Also what type of
> contact used in Yade? I think the same like in Jan's article. May be you
> have something else Jan, especially for my type of simulation?

Well, the described contact law is "standard", one of the most basic linear
elastic laws. The same approach is used in CpmMat. The "strain" I describe
in the article is strain of interaction, relative displacement and rotation
of particles from the point of view of the interaction, completely
different strain than that you want to plot :-)

How to plot contact forces and how much vectors of force will have each
> sphere?

there is VTKExporter.exportInteractions. Or you can use VTKRecorder with
recorders=['intr']

Jan, I'm not happy yet with middle stresses because they are bigger:) And
> by the way what do you think about Bruno's proposal to devide macroStress
> = sphereStress*sphereVolume/cubeVolume and what should i do in the
> non-regular packing?
> As said in the future there will be non-regular areas inside.

you have to define sphere's volume in a different manner then just cube
with edge 2*r :-) e.g. volume of respective Voronoi cell

I'v already read Yade manual with DEM background but there is no simple
> formulas for strain and stress, for example in you article Jan strain
> described much better. Also it's interesting to read about interactions
> laws in CpmMat material model which i used now.

Again, in my article I described idfferent strain :-) there is explanation
in User's manual [3] with a reference (the article is available on
Reseachgate), equations 39 and 40.

Contact law is desribed in [4] with
Kn = young*A/L, A=pi*r*r for unisize particles, L being particles' centers
distance
Kt = Kn*poisson

So to conclude the main things:
- please ask in another questions concentrating on one specific problem
- the value of forces you applied on particles now are not correct (does
not correspond to defined pressure)
- bodyStressTensors is ok for inner particles, not so ok for boundary
particles without prescribed force and definitely not ok for particles with
prescribed force, you can compute their stress "manually"
- using smaller and smaller particles and the same parameters should
converge to one value (with correctly applied forces)

cheers
Jan

[1]
https://yade-dem.org/doc/yade.export.html#yade.export.VTKExporter.exportInteractions
[2]
https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.VTKRecorder.recorders
[3] https://yade-dem.org/doc/user.html#micro-stress-and-micro-strain
[4] https://yade-dem.org/doc/formulation.html#stress-evaluation-example

2015-07-16 22:36 GMT+02:00 Alexander <email address hidden>:

> Question #269063 on Yade changed:
> https://answers.launchpad.net/yade/+question/269063
>
> Alexander gave more information on the question:
> Also the last MWE:
>
>
> import __builtin__
> from yade import export
>
> ###################################################
> # define materials and model configuration
>
> E = 5.2e11#0.65e11#1.338e11#2e11 # Young's modulus of model need to be
> calibrated
> v = 0.3# Poisson's ratio
> P = 150e6 # initial pressure value
> d = 7850 # density
>
> r = 0.25 # spheres radius
>
> # Enlarge interaction radius between spheres using "interaction_radius"
> parameter (for example in uniax.py this value is 1.5)
> interaction_radius = 1.5
>
> F = P*r
>
> # define plate material, create "dense" packing by setting friction to
> zero initially
> O.materials.append(CpmMat(young=E,
> frictionAngle=0,
> poisson=v,
> density=d,
> sigmaT=3.5e6,
> epsCrackOnset=1e-4,
> neverDamage=True,
> isoPrestress=0,
> relDuctility=30,
> label = 'mat'))
>
> # represent plate like a set of regular monosized set of spheres
> # also set boundary conditions via predefined tensile force for spheres on
> ABCD and
> # fixed spheres on KGHO
>
> O.bodies.append(pack.regularOrtho(pack.inAlignedBox((0,0,0),(16.1,16.1,2.1)),r,0))
> for b in O.bodies:
> id = b.id
> if b.state.pos[1]+r > 15.9:
> O.forces.addF(id,(0,F,0),permanent=True)
> if b.state.pos[1]-r < 0.1:
> O.forces.addF(id,(0,-F,0),permanent=True)
>
>
> ###################################################
> # define engines
>
> # simulation loop
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interaction_radius,label='bo1s')]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interaction_radius,label='ig2ss')],
> [Ip2_CpmMat_CpmMat_CpmPhys()],
> [Law2_ScGeom_CpmPhys_Cpm()]
> ),
> CpmStateUpdater(realPeriod=1),
> NewtonIntegrator(damping=0.4)
> ]
>
> ###################################################
> # start simulation and compute strain and stress
>
> # try to run script with qt graphical interface
> try:
> yade.qt.Controller(), yade.qt.View()
> except:
> print 'Qt graphical interface is not avaliable'
>
> # set the integration timestep to be 1/2 of the "critical" timestep
> O.dt=.5*utils.PWaveTimeStep()
>
> # run one single step
> O.step()
>
> # reset interaction radius to the default value
> bo1s.aabbEnlargeFactor=1.0
> ig2ss.interactionDetectionFactor=1.0
>
> # run simulation, until static equilibrium will not reached
> while unbalancedForce()>1e-2:
> O.run(10,True)
>
> # compute stress tensor for each body
> stresses = bodyStressTensors()
>
> ###################################################
> # save data to vtk. file
> __builtin__.my_stresses = stresses
> vtk = export.VTKExporter('result')
>
> vtk.exportSpheres(what=[('radius','b.shape.radius'),('displacement','b.state.displ()'),('stress','my_stresses[
> b.id]')])
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

>So is it correct to apply force (pressure) to sphere without taking into account it's size, or i need to multiply by it's radius, diameter or volume?

None of them ;)
You need to multiply pressure by area, as always. When you change the size of spheres you also change the thickness of the plate.

>Also Bruno about your formula macroStress = sphereStress*sphereVolume/cubeVolume. How it will work for non regular packing?

https://yade-dem.org/doc/user.html?highlight=microstress#micro-stress-and-micro-strain

>Yes Ansys solution is exact for uniaxial tension

I don't think Ansys gives the displacement field that you can derive analyticaly for a sphere array. It gives the exact solution for a homogeneous isotropic material, instead. I really mean to write on the paper the displacement of each sphere in an array, you don't need DEM for that.
Exemple for a single chain of spheres (diameter D) where sphere 0 is fixed, the displacement of sphere "i" is:
u[i]=F/kn*i

An equivalent expression exhibits equivalent stress and strain:
u[i]/(iD)=F/(kn*D) = F/D² * 1/Young
since kn (contact stiffness) is defined as kn=Young*D.
The meaning of this last one is strain=stress/Young

You should go through this type of thing to understand how DEM will or will not do what you expect. It will also help clarify what is strain and what is stress.

Besides I vote as Jan suggests about closing this discussion.

Hello Jan and Bruno.

I'm so sorry about mess with the applied pressure:( Of course i gotta use force F=P*2*r*2*r for each sphere where P is predefined pressure and r - is sphere radius. The second variant is define total_force = P * surface_area and then divide total_force by the number of boundary particles where force is applied.

So as you said i close this discussion now. But don't forget about me, i'll reopen it with the problem of stress and strain computation with the name "Stress/strain evaluation in the uniaxial tension test":)

Excuse me but i wanna ask one question before open new discussion:

>- bodyStressTensors is ok for inner particles, not so ok for boundary
>particles without prescribed force and definitely not ok for particles with
>prescribed force, you can compute their stress "manually"

Could i just create some additional boundary layer of spheres surrounded the plate and then delete it after simulation has been done? So in this case spheres with prescribed force and boundary spheres will be deleted and only inner spheres will remain. So for them bodyStressTensors is ok.

 Thank's a lot for given advise, especially Jan , but Bruno also hepled well:)

with best regards, Alexander

>
>
> Could i just create some additional boundary layer of spheres
> surrounded the plate and then delete it after simulation has been done?
> So in this case spheres with prescribed force and boundary spheres will
> be deleted and only inner spheres will remain. So for them
> bodyStressTensors is ok.
>
>
sure :-) e.g. VTKExporter offer to specify only some spheres to export.

But to be scietifically correct, I would let the results as they are,
commenting why the boundary is different and tell that in limit this
anomally is a set of measure zero :-) Or you can postprocess the results a
bit (e.g. taking the prescribed force into account for stress computation)

cheers
Jan

Thank's a lot again

So close this discussion.

with regards, Alexander